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Abstract 

Precise  position  determination  and  control  is  necessary  to  accomplish  proposed 
satellite  formation  flying  missions  of  ground  movement  target  indication  and  synthetic 
aperture  radar.  This  thesis  combines  the  estimation  and  control  techniques  of  past  AFIT 
theses  with  various  time-varying  and  time-invariant  LQG  control  methods.  Linear  time- 
invariant  control  is  ideal  for  on-board  satellite  estimation  and  control  applications, 
freeing-up  the  satellite’s  limited  computational  capacity.  Using  a  dynamics  frame 
transformation  from  the  nodal  frame  to  an  orbital  frame,  a  higher  fidelity,  time-periodic 
model  produced  nearly  identical  results  for  either  time-varying  or  time-invariant  control 
for  many  scenarios.  Scenarios  included  initial  perturbations  in  the  radial,  in-track,  and 
cross-track  directions  as  well  as  increased  magnitude  perturbations;  step  size  increase 
from  0.2  seconds  to  2  seconds;  and  increased  and  reduced  measurement  noise  level 
scenarios  versus  the  standard  absolute  GPS  receiver  noise  level.  Results  obtained 
indicate  the  ability  to  control  within  the  error  range  of  the  measurements  (centimeter- 
level  and  better)  using  realistic  noise  and  dynamics  perturbations. 
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PERIODIC  METHODS  FOR  CONTROLLING  A  SATELLITE  IN  FORMATION 


I.  Introduction 


1.1  Background 


In  tackling  this  subject,  two  questions  come  to  the  forefront.  Why  focus  on  time 
periodic  methods  and  why  satellite  formations? 

Periodic  characteristics  are  observed  both  naturally  and  artificially  in  many 
mechanisms,  most  notably  in  satellite  orbits  and  orbit  perturbations.  Since  orbital  control 
usage  is  a  big  driver  in  the  lifetime  of  all  satellites,  only  using  control  when  necessary 
should  be  the  goal  of  the  control  system.  If  perturbations  are  periodic  and  small  enough 
that  they  can  be  ignored,  there  is  no  reason  to  spend  fuel  to  control  them.  The  periodic 
nature  of  the  perturbation  will  correct  itself.  The  detail  of  how  small  the  perturbation  can 
be  before  it  warrants  correction  must  be  defined  and  justified  via  simulations  before  it 
can  be  used  for  actual  flight  tests. 

The  importance  of  modeling  the  periodic  characteristics  is  multiplied  as  future 
space  missions  will  use  several  small  satellites  flying  in  formation  to  serve  many  of  the 
same  purposes  as  a  single  large  satellite. 


Large,  specialized  satellites 


Paradigm  Shift 


Virtual  satellite  - 

Clusters  of  cooperating  satellites 
Figure  1.  Paradigm  Shift  [1] 
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This  paradigm  shift  will  produce  many  rewards:  lower  overall  mission  cost, 
increased  system  reliability,  improved  science  return,  extended  mission  life,  and 
expandability.  Satellite  formation  has  been  a  widely  studied  area  for  nearly  a  decade, 
with  a  patent  for  the  concept  granted  in  1983  [24].  Many  space  agencies  currently  have 
actual  or  proposed  satellite  formations  on  their  drawing  boards  or  operating  in  space. 

The  European  Space  Agency’s  Cluster  II  project,  launched  in  2000,  uses  four 
satellites  in  tetrahedral  formation  to  study  electrical  and  magnetic  phenomena  inside  the 
magnetosphere  and  in  the  solar  wind  upstream  of  the  Earth. 


Figure  2.  ESA’s  Cluster  II  and  NASA’s  StarLight  Satellites 


NASA’s  StarLight  mission,  consisting  of  two  satellites  in  formation  to  be 
launched  in  2005,  will  be  the  first  spacebome  stellar  interferometer.  By  flying  in 
formation,  Starlight’s  two  small  telescopes  will  achieve  the  angular  resolution  (in  one 
direction)  of  a  regular  telescope  mirror  125  meters  in  diameter. 

The  United  States  Air  Force  (USAF)  has  focused  its  efforts  on  the  Technology 
Satellite  of  the  2 1st  Century  (TechSat  21).  The  Air  Force  Research  Laboratory  (AFRL) 
leads  the  project. 
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Figure  3.  AFRL’s  TechSat  21  Satellite  [1] 


As  AFRL’s  website  explains,  the  initial  focus  of  TechSat  21  is  a  space-based 
radar  mission  for  Ground  Moving  Target  Indication  (GMTI)  employing  a  formation  of 
eight  micro-satellites  [1].  This  program  hopes  to  achieve  three  times  the  capability  of 
current  assets  at  one-third  their  cost. 
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Table  1.  AFRL’s  TechSat21  Cost  Savings  [T 


Conventional 

Evolutionary 

Revolutionary 

Rigid  Panel  Deployable 

Integrated  Antenna 

Distributed 

Spacecraft 

UH£> 

mg 

— 

Satellite  Mass 

12,000  kg 

4,400  kg 

Eight  100  kg 
satellites  per  cluster 

Age  of  Technology 

1 996  Technology 

2003  Technology 

2005  Technology 

Number  of  Satellites 
Required 

40  satellites 

40  satellites 

40  satellite  clusters 

Launch  Vehicle 

Titan  IV 

Atlas  II 

Taurus/Athena  II 

Launch  Cost 

$300M 

$90M 

$25M 

Total  System  Cost 
(10  year  life  cycle, 
normalized) 

1.0 

0.6 

0.3 

The  reduced  life  cycle  cost  will  be  achieved  by  reduced  launch  costs  and  overall 
manufacturing  cost  through  mass-produced  identical  micro-satellites.  The  TechSat  21 
program  manager  even  speculates  the  possibility  that  TechSat  21  satellites  could  be 
launched  from  operational  military  jets  such  as  the  F-15  [23]. 

The  GMTI  mission  of  detecting  ground  targets  in  a  military  theater  for  the  USAF 
is  currently  being  accomplished  by  Joint  Surveillance  and  Targeting  Attack  Radar 
System  (JSTARS)  aircraft.  The  proposed  Techsat  21  system  could  move  this  GMTI 
mission  from  an  airborne  platfonn  to  space.  Better  performance  will  come  from 
extremely  large  effective  aperture  sizes  and  multi-mission  capability.  For  radars,  large 
apertures  improve  performance  with  better  slow  target  detection  and  better  resolution.  In 
fact  if  the  range  to  the  target,  R,  and  the  radar’s  wavelength,  X  ,  are  assumed  constant,  the 
resolution  of  an  object  on  the  ground,  r,  is  inversely  proportional  to  S,  the  synthetic 
aperture  size  [16]. 


R  ■  A, 

r  = - 

2  ■  S 


(1) 
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The  key  to  the  concept  is  the  close  formation  of  microsatellites,  each  with  a 
receiver  that  detects  both  the  return  from  its  own  transmitter  and  the  bistatic  responses 
from  the  orthogonal  transmit  signals  of  the  other  satellites  in  its  cluster.  Pulses  are 
separated,  and  range,  angle  of  arrival,  and  Doppler  processing  is  computed.  The 
collaborating  small  apertures  of  TechSat  21  can  act  as  one  large  array  which  can  be 
reconfigured  to  satisfy  and  optimize  requirements  for  a  variety  of  missions. 

Passive  Radiometry  Mission 


(Narrow  beam/wide  area  coverage) 

Figure  4.  AFRL’s  TechSat  21  Missions  [1] 

Possible  autonomous  capability  would  allow  the  satellite  cluster,  working  in 
conjunction  with  the  Global  Positioning  System,  to  automatically  reshape  its  formation  to 
optimize  its  geometry  based  on  the  mission  assigned  to  the  cluster  by  the  ground 
commander.  For  example,  a  “geo-location”  mission  such  as  finding  a  downed  pilot 
would  require  that  one  satellite  in  the  cluster  move  out  from  the  group  to  a  distance  of 
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approximately  eight  kilometers.  Or  for  spotting  moving  ground  targets,  closer  intervals 
on  the  order  of  several  hundred  meters  would  be  required  [23]. 

Improved  reliability  will  be  achieved  with  graceful  degradation  of  capability  and 
reconfigurability  to  minimize  the  effects  of  a  single  failure.  The  satellite  clusters  also 
allow  an  inherent  adaptability  easing  performance  upgrades  by  allowing  an  upgraded 
satellite  to  join  a  cluster  vice  replacing  another  satellite.  Much  like  a  network  of  desktop 
computers  linked  together,  they  can  be  expanded  as  needed  without  throwing  out  the 
entire  system. 

“The  availability  of  such  low  cost,  highly  capable  microsatellites  will  lead  to  new 
ways  of  doing  business  in  space,  and  this  ‘personal  satellite’  paradigm  can  revolutionize 
the  space  industry  in  much  the  same  way  as  the  personal  computer  has  done  to  the 
computer  industry,”  said  Dr.  Alok  Das,  chief  technologist  for  AFRL’s  Space  Vehicles 
Directorate  [12]. 


1.2  Past  Research  -  Standing  on  the  Shoulders  of  Giants 

1.2.1  Satellite  Cluster  Research  at  AFIT  from  1988  -  1991.  A  number  of 
AFIT  theses  in  the  past  have  utilized  Kalman  filter  techniques  to  study  relative  position 
errors  within  a  satellite  cluster  [19,  24],  The  dynamics  were  based  on  Clohessy-Wiltshire 
equations  with  near  circular  orbits.  Ward  used  a  point  mass  orbit  in  his  truth  model,  but 
Middendorf  included  a  more  realistic  truth  model.  He  included  perturbations  in  the  J2 
tenn  of  the  Earth’s  geopotential.  For  relative  positioning  they  both  used  a  U-D 
covariance  factorization  version  of  the  Kalman  filter  comparing  range  data  between 
satellites  synchronized  with  clock  pulses.  They  achieved  a  relative  position  error  of  3 
cm.  This  relative  centimeter-level  accuracy  was  achieved  without  the  use  of  Global 
Positioning  System  (GPS)  satellite  positioning  data. 
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1.2.2  Satellite  Cluster  Research  at  AFIT  from  2000  -  2001.  More  recent  AFIT 


theses  incorporated  GPS  positioning  data  into  a  Bayes  filter.  Davis  attempted  to  extend 
from  relative  position  determination  of  the  past  to  include  absolute  position 
determination,  but  he  encountered  numerical  difficulties  due  to  the  extreme  accuracy 
requirements  of  his  filter  [9].  The  double  precision  computational  capability  was  not  able 
to  accurately  generate  covariance  inversions  for  his  filter.  Bordner  revisited  the  previous 
theses  and  focused  mainly  on  relative  positioning  [4].  Following  a  new  approach  by 
Wiesel  [28],  Bordner  updated  the  dynamics  model  by  solving  the  time  periodic, 
linearized  system  using  Floquet  Theory.  More  accurate  than  the  Clohessy-Wiltshire 
solutions  used  in  previous  theses,  this  dynamics  model  included  all  zonal  harmonics  of 
the  Earth  to  include  sectoral,  tesseral,  and  air  drag  perturbations.  Bordner’s  dynamics 
truth  model  was  also  used  in  this  current  thesis  as  discussed  in  Chapter  III.  Bordner 
additionally  expanded  on  the  method  of  relative  data  collection  from  previous  theses.  He 
compared  results  using  the  relative  range  information  from  satellite  timing  signals,  as 
was  the  basis  for  the  AFIT  theses  from  a  decade  ago,  with  both  absolute  GPS  data  and 
Carrier  Phase  Differential  GPS  data.  Using  Differential  GPS  data,  Bordner  achieved  a 
relative  position  error  of  2  cm. 

Another  AFIT  thesis  authored  by  Irvin  studied  several  linear  and  nonlinear 
feedback  control  methods  for  satellite  formation  reconfigurations  [9].  One  of  his  more 
successful  methods,  the  state  dependent  Riccati  equation,  was  used  as  a  starting  point  for 
the  regulator  in  the  second  phase  of  this  thesis.  This  current  thesis  diverged  from  Irvin’s 
work  by  using  the  assumed  superior  dynamics  model  as  used  by  Bordner  in  place  of 
Irvin’s  Clohessy-Wiltshire  method.  The  combined  filter  and  regulator  for  this  thesis 
should  provide  a  good  combination  building  on  Bordner’s  estimation  work  and  Irvin’s 
control  method. 
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1.3  Objectives 

In  order  for  the  TechSat  21  project  to  work  and  be  cost  effective,  AFRL  lists  a 
number  of  fundamental  questions  that  must  be  answered  [1].  The  focus  of  this  thesis  is  to 
answer  “what  orbits  or  orbital  control  techniques  are  possible  to  maintain  the  desired 
configuration?” 

This  thesis  will  look  at  two  different  periodic  frames  of  reference  for  the 
dynamics  of  the  system,  the  nodal  frame  and  the  orbital  frame.  In  each  of  these  frames, 
time-varying  and  time-invariant  control  techniques  will  be  compared.  Variations  on  the 
system  model  will  include  scenarios  with:  initial  deflections  in  each  direction,  increased 
magnitude  of  initial  deflection,  reduced  step  size,  decreased  measurement  noise,  and 
increased  measurement  noise.  Some  of  these  scenarios  may  improve  performance  of  the 
system,  but  taken  as  a  whole,  all  scenarios  will  test  the  fidelity  of  each  model  for  many 
environments. 

The  ideal  result  of  this  thesis  would  be  to  create  a  robust  output  feedback 
compensator  that  handles  a  wide  variety  of  realistic  perturbation  scenarios  minimizing 
both  fuel  usage  and  Root  Mean  Squared  (RMS)  error.  The  stringent  space-based  radar 
application  would  require  position  accuracy  to  one  tenth  of  the  radar’s  wavelength  [8], 
thus  requiring  millimeter-level  accuracy.  While  past  theses  have  shown  that  positions 
can  be  estimated  to  the  centimeter  level,  we  would  also  like  to  be  able  to  control  our 
satellite  to  the  centimeter  level  or  better.  This  may  not  be  an  easy  task  given  the 
preliminary  assumed  conditions  of  estimation  using  absolute  GPS  position  data 
constrained  with  over  one  meter  uncertainty. 


II.  Dynamics  Theory  and  Methodology 


2.1  Satellite  Relative  Motion 

As  mentioned  in  the  previous  chapter,  this  current  thesis  used  the  same  method  as 
Bordner  for  calculating  the  dynamics  for  the  truth  model  [4].  Relative  satellite  motion 
has  often  been  studied  and  applied  using  satellite  rendezvous  dynamics  models  based  on 
Clohessy-Wiltshire  (CW)  equations  with  near  circular  orbits.  Clohessy  and  Wiltshire 
followed  on  from  the  work  of  Hill  from  the  1800’s.  They  reduced  the  dynamics  of 
satellite  relative  motion  into  a  two-body  problem.  They  assumed  a  circular  orbit, 
linearized  the  system  of  equations,  and  produced  a  solution.  Although  this  solution  was 
originally  intended  to  describe  the  relative  motion  of  a  closing  satellite  docking  with  a 
target  satellite  or  space  station,  CW  equations  are  usually  chosen  as  the  place  to  start  for 
satellite  formation  analysis. 

2.1.1  Clohessy-Wiltshire.  While  the  CW  method  was  not  used  for  this  thesis,  it 
was  useful  to  look  at  its  setup  and  compare  it  with  the  method  employed  by  Bordner  in 
his  thesis  [4].  The  CW  equations  in  their  simplest  form  did  not  model  any  perturbations 
in  their  solution  to  the  linearized  circular  orbit.  For  our  earth-orbiting  satellites,  the  CW 
equations  are  expressed  in  differential  form 

x  -2  ■  n  ■  y-3  ■  n2  •  x  =  0 

y  +  2  ■  n  •  x  =  0  (2) 

z  +  n2  ■  z  =  0 

where  n  is  the  mean  motion 
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n  = 


(3) 


3  2* 

where  jj.e=3 .986E 14  m  /s“  is  the  earth  gravitation  parameter  and  Rsat  is  the  radius 
from  the  center  of  the  earth  to  the  satellite  in  orbit.  The  standard  Clohessy-Wiltshire 
reference  frame  employs  x,  y,  and  z  corresponding  to  radial,  in  track,  and  orbit  normal. 
For  the  “free  motion”  CW  model,  Equations  2  were  all  set  equal  to  zero.  More  advanced 
CW  models  would  include  perturbations  added  on  the  right  hand  sides  of  these  equations. 


Figure  5.  Clohessy-Wiltshire  Reference  Frame 


Following  a  free  motion  CW  model,  satellites  would  be  constrained  to  an  ellipse 
of  2  to  1  dimensionality  along  the  trajectory  of  the  orbit  of  the  satellite  and  in  the  vertical 
plane  of  motion  [21].  The  actual  dynamics  of  the  system  would  be  periodic,  and  even 
though  the  relative  distances  between  satellites  would  remain  variable,  the  overall  size  of 
the  cluster  would  be  bounded.  Sedwick  further  described  the  destruction  of  the  free 
motion  CW  model  due  to  the  overwhelming  nature  of  the  real  world  secular  and  periodic 
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perturbations.  The  goal  for  satellite  formation  modeling  would  be  to  include  all  forces 
that  cause  the  orbital  position  to  diverge.  Benign  perturbing  forces  driving  satellites  to 
move  periodically  or  in  unison  could  remain  un-modeled. 

But  even  adding  significant  perturbation  estimates  to  the  CW  equations  would  not 
be  expected  to  supply  the  best  estimates  for  satellite  fonnation  navigation  [4].  Hill’s 
equations  were  based  on  a  time-invariant  system,  which  limits  the  credibility  of  the 
dynamics  for  its  application  for  periodic  systems  such  as  satellite  formations.  The  CW 
equations  were  designed  for  satellite  rendezvous,  not  satellite  formations.  A  new  relative 
satellite  dynamics  model  was  needed  to  adequately  cover  the  nature  of  the  satellite 
formations. 

2.1.2  Wiesel.  As  mentioned  above,  this  current  thesis  used  the  same  method  as 
Bordner  for  calculating  the  dynamics  for  the  truth  model  [4].  This  method  was  a  result  of 
work  accomplished  by  Dr.  William  Wiesel  in  the  area  of  relative  satellite  motion 
developed  to  meet  the  severe  relative  positioning  accuracy  requirements  inherent  in 
space-based  radar  applications  [27,  28].  His  method,  similar  to  the  classic  Clohessy- 
Wiltshire  solution  for  relative  satellite  motion,  took  a  new  approach  by  linearizing  about 
a  nearly  circular  periodic  orbit  as  a  reference  solution.  The  reason  for  integrating  a 
periodic  orbit  versus  some  other  particular  orbit  was  that  the  orbital  behavior  was  better 
understood  as  time  drove  towards  infinity  [26].  This  was  because  in  the  inertial  frame 
the  orbit  behavior,  given  the  same  initial  conditions  and  perturbations,  would  repeat 
every  period  plus  an  orbital  regression  along  the  equator.  To  work  with  an  orbit  that 
closed  on  itself  exactly,  the  nodal  frame  was  introduced.  This  allowed  a  resulting  linear, 
time  periodic  system  instead  of  a  linear,  time  invariant  system. 

WieseTs  process  followed  a  Hamiltonian  method  [4,  28],  A  Hamiltonian  function 
was  used  to  describe  the  orbital  motion  of  a  small  body  traveling  around  the  earth.  With 
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dimensionless  units,  setting  the  radius  of  earth  Re  =  1,  gravitational  constant  G  =  1,  and 
the  mass  of  earth  Me  =  1 ,  the  inertial  frame  Hamiltonian  function  became 


The  inertial  state  vector  was  defined  as  IT  =  (X,Y,Z,Px,Py,Pz),  and  the  radius  r  was 
defined  as  r  =  ylx2  +  Y2  +  Z2  .  The  P;  momenta  states  were  defined  as  inertial  velocity 
components  because  the  Hamiltonian  is  in  terms  of  mass  per  unit  satellite.  An  advantage 
of  using  the  Hamiltonian  method  was  the  convenience  of  adding  perturbations. 
Additional  conservative  perturbations  could  be  modeled  by  including  the  potential 
function  of  the  perturbing  force.  Non-potential  forces  such  as  air  drag  were  expanded 
about  the  periodic  orbit,  and  then  added  to  the  Hamiltonian  equations  of  motion  for  the 
momenta  states. 

In  the  inertial  frame,  the  actual  orbit  would  regress  along  the  equator  and  not 
return  to  the  same  position  after  one  period  [28],  Because  of  this  regression,  the  system 
was  transfonned  to  a  nodal  frame  where  the  orbit  would  actually  close  on  itself  thus 
becoming  periodic.  The  resulting  Hamiltonian  became 

h' = \  (p;1 + p; + p;= )+ (p;  •  v  -p;  •  x  )+ vfz'.o  o 

where  the  prime  superscript  symbol  denoted  the  state  vector  components 
transformed  into  the  nodal  frame,  N.  The  component  v(z  , r)  included  the  conservative 

zonal  hannonics  terms  transfonned  into  the  nodal  frame.  The  nodal  frame  state  vector 
was  defined  NT  =  ( X  ,  Y  ,  Z  ,  P’ ,  P| ,  P’).  The  beauty  of  this  approach  was  that  the  nodal 

frame  was  not  the  only  frame  in  which  the  system  was  periodic.  The  nodal  system  could 
be  moved  back  to  the  classic  CW  frame  (Figure  5)  and  still  keep  its  periodic  quality.  To 
accomplish  this,  the  periodic  orbit  was  initially  transfonned  to  the  origin  of  the  current 
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coordinates.  The  periodic  orbit  was  subtracted  off  from  the  nodal  state  vector,  giving 
new  coordinates 

x  =  X'-X'0 


y  =  Y'-Y0' 


z  =  z'-z; 

Px  =Pxo(t) 


(6) 


Py 


Pz 


=p;0(t) 


where  the  subscript  “o”  denoted  the  periodic  orbit.  This  transfonnation  produced 
the  Hamiltonian  with  its  second  order  and  above  terms  evaluated  on  the  periodic  orbit. 
Using  the  new  coordinate  vector  LT  =  (x,y,z,px,py,pz)  defined  by  Equation  6,  the 
Hamiltonian  now  became 

H'  =  E(Hi-L„-L|i)+i.(H^.L„.Lf|.LT)+-  (7) 

The  H^  and  H''Py  components  were  fully  symmetric  partial  derivative  tensors 
computed  for  the  periodic  orbit  defined  as 


32H' 


H 


“l”  3N„3Nf3N, 


(8) 


where  repeated  Greek  indices  were  a  summation  from  one  to  six. 
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A  final  transformation  rotated  the  solution  into  the  classic  CW  frame  (radial,  in¬ 
track,  and  cross-track).  To  generate  a  linear  system,  the  Hamiltonian  K  was  truncated  to 
include  only  terms  following  the  second  order 

K;  =  (h;,  R$  ■  Rg> )+  P  •  (r  •  RT  )•  q  <9> 

where  H*fl  was  defined  above  in  Equation  8,  R  ^  and  Rj^  were  6  by  6 


dimensioned  block  diagonal  matrices  with  two  copies  of  the  nodal  to  orbital  frame 
rotation  matrix  R  on  the  diagonal,  the  Z  components  were  the  canonical  state  vectors  in 
the  orbital  frame,  and  the  p  and  q  components  were  3  by  1  vectors  holding  the  momenta 
and  coordinates  from  the  local  frame  state  vector  L  [28],  The  canonical  equations  of 
motion  would  then  be  calculated  by  using  Hamilton’s  Equations 


q, 


P,  = 


3K2 

aPl 

3q, 


(10) 


to  evaluate  the  Hamiltonian  function  in  Equation  4. 


2.1.3  Floquet  Theory.  The  equations  of  motion  calculated  in  the  previous 
section  can  be  integrated  for  one  period  to  obtain  the  desired  periodic  orbit  [4,  28], 
Linearized  equations  are  created  by  expanding  about  this  periodic  orbit.  The  periodic 
nature  of  these  equations  made  the  approach  to  solve  them  a  bit  more  complicated  than 
the  simple  CW  solution.  The  solution  was  based  on  Floquet  theory  that  linearizing  about 
a  periodic  orbit  resulted  in  a  periodic  coefficient  linear  system.  The  time  periodic  linear 
differential  equations  were  of  the  form 

dx  =  A(t)-3x  (11) 

Floquet  found  that  the  solution  to  Equation  1 1  was  the  state  transition  matrix 
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<J>(t,0)  =  F(t)  ■  eJt  ■  F_1  (0) 


(12) 


where  F(t)  is  a  periodic  matrix,  and  J  is  a  matrix  of  system  frequencies  in  the 
Jordon  normal  form.  The  J  matrix  is  typically  diagonal  with  Poincare  exponents  as 
elements  of  the  diagonal.  This  Floquet  state  transition  matrix  from  Equation  12  is  similar 
to  the  constant  coefficient  form,  but  it  is  unique  in  its  periodic  nature.  To  find  the 
solution  for  <f>(t,0) ,  both  J  and  the  periodic  matrix  F(t)  needed  to  be  calculated.  Equation 

1 1  was  integrated  for  one  period  to  calculate  the  mondromy  matrix.  For  a  periodic  F(t), 
F(t  +  x  )  =  F(t).  So  plugging  in  t  =  x  ,  Equation  12  became 

<J>(x,0)  =  F(x)  ■  eJx  •  F_1  (0)  (13) 

after  rearranging  terms,  F(0)  was  recognized  as  the  eigenvector  matrix  for  <t>(x,0) 

F(0)  ■  <f>(x,0)  ■  F_1  (0)  =  ejT  (14) 

The  eigenvalues  of  <f>(x,0)  ,  were  calculated  at  F(0)  using  the  relation 

A.;  =  e“‘  ■  x  (15) 

where  co,  were  the  Poincare  exponents  of  the  J  matrix.  Rearranging  this  equation 

to  solve  for  CO;  would  eventually  lead  to  the  solution  for  J. 

CO;  =  — lnA.;  (16) 

x 

Substituting  Equation  12  into  Equation  1 1 

F  =  A(t)  •  F(t)  -  F(t)  ■  J  (17) 

Integration  of  Equation  17  for  one  period  provided  a  solution  for  <J>(t,0) .  The  end 

result  was  an  inertial  position  vector  for  the  satellite 

I(t)  =  R<2)  ■  (N0 (t)  +  (R(2>)t  ■  F(t) ■  eJt  •  F-‘(t0) •  Z(t0 ))  (18) 

and  the  satellite’s  position  relative  to  the  periodic  orbit  in  the  orbital  frame  was 
Z(t)  =  F(t)  ■  eJt  ■  F“'(to)  •  Z(t0)  (19) 
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2.1.4  Wiesel  vs.  Clohessy-Wiltshire.  Wiesel’s  dynamics  model  approach 
attempted  to  better  predict  satellite  behavior  and  minimize  the  fuel  required  to  control  the 
satellite  by  modeling  all  major  perturbations.  As  Bordner  showed  in  his  thesis,  this 
method  improved  on  the  Clohessy-Wiltshire  approach  in  three  areas  [4],  First,  Wiesel’s 
method  inherently  retained  more  satellite  behavior  information  over  time,  thus  it 
introduced  less  error  into  the  dynamics  model.  Second,  two  secular  modes  (rather  than 
the  usual  one  secular  mode)  were  modeled  by  the  dynamics.  These  secular  modes  would 
be  better  compensated  for  without  using  extraneous  fuel  to  try  to  negate  them.  Third, 
Wiesel  included  a  better  representation  of  the  earth’s  zonal  harmonics  (zonal  effects  were 
modeled  to  the  14th  order  in  this  thesis),  air  drag,  and  gravitational  harmonics.  Better 
modeling  of  perturbations  produced  a  more  realistic  model  to  build  upon.  This  solution 
improved  on  the  Clohessy-Wiltshire  solution  by  three  orders  of  magnitude  [28]. 

2.2  Actual  Thesis  Dynamics  Setup 

The  first  step  in  setting  up  the  estimator/controller  was  to  define  where  you  are. 
As  mentioned  earlier  in  this  chapter,  the  standard  Clohessy-Wiltshire  reference  frame 
was  used  with  x,  y,  and  z  corresponding  to  radial,  in  track,  and  orbit  normal  as  is  shown 
in  Figure  5.  This  earth-centered  inertial  reference  frame  was  defined  in  terms  of  distance 
units  (DU),  with  one  DU  approximately  equal  to  the  equatorial  radius  of  the  earth, 
6378.145  km,  and  time  units  (TU),  with  one  TU  approximately  equal  to  13.447  minutes 
[2].  Based  roughly  on  Techsat  2 1  planned  characteristics,  the  orbit  height  was  set  at  0. 1 
DU  (637.8  km)  with  a  7.244  TU  (97  minute)  period  for  a  near  circular  orbit.  An  existing 
program  was  used  to  generate  the  initial  velocities  corresponding  to  the  chosen  orbit. 

The  initial  state  was  defined  as  a  vector  of  positions  and  velocities: 
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(20) 


X 

l.l(DU) 

y 

O.O(DU) 

z 

O.O(DU) 

X 

-0.0012936533(DU/TU) 

y 

0.5 1 53407938(DU  /  TU) 

z 

0.802595733 1(DU/TU) 

2.2.1  Orbit  Propagation.  Following  the  method  Bordner  used  for  his  thesis,  the 
state  was  propagated  from  initial  conditions  through  one  orbit  using  a  program  written  in 
the  fortran  computer  language.  The  source  code  is  primarily  written  by  Dr  Weisel  with 
some  additions  to  meet  the  specific  requirements  of  this  thesis.  The  main  program, 
Po.f90,  uses  two  main  subroutines,  haming  and  rhs.  Subroutine  haming  is  an  ordinary 
differential  equations  integrator  using  a  fourth  order  predictor-corrector  algorithm. 
Haming  collects  the  last  four  values  for  the  state  vector  and  extrapolates  them  to  predict 
the  next  value.  It  then  corrects  the  extrapolated  prediction  to  find  a  new  value  for  the 
state  vector.  Subroutine  rhs  calculates  the  equations  of  motion  calling  on  the  EGM96 
earth  gravity  model  and  geopoten  subroutines  to  add  to  the  system  a  more  accurate  model 
for  the  zonal  harmonics  as  well  as  air  drag  and  gravitational  harmonics. 

It  is  worth  noting  that  an  intrinsic  part  of  this  thesis  is  the  periodic  nature  of  the 
orbit.  Once  the  state  vector  is  propagated  over  one  orbit,  the  state  vector  (and  its 
associated  A  matrix)  will  repeat  exactly  over  following  orbits.  The  word  “exactly”  is 
used  in  a  perfect-world  sense.  The  perturbing  effects  that  cause  the  dynamics  not  to 
repeat  “exactly”  will  be  modeled  as  either  dynamics  noise  (w)  or  measurement  noise  (v) 
to  the  system.  For  this  thesis,  the  orbit  closed  on  itself  (initial  state  equals  final  state) 
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with  errors  on  the  order  of  IE- 1 1  DU  and  DU/TU.  This  error  was  comparable  to  about 


IE-5  meters  and  IE-8  meters/second. 


2.3  Frame  Variations  for  the  A  matrix 


Along  with  the  state  vector  propagated  over  the  orbit,  the  fortran  program  also 
propagated  the  A  matrix  associated  with  the  state  vector.  This  A  matrix  was  derived  as 
the  time  periodic  solution  to  the  uncontrolled  dynamics  problem 

X(t)  =  A(t)  -X(t)  (21) 

The  original  method  of  computing  the  A  matrix  was  in  the  nodal  frame  as 
described  in  a  previous  section.  Within  an  orbital  period,  elements  of  this  A  matrix 
would  vary  up  to  an  order  of  magnitude  1E1.  Using  a  transfonnation  supplied  by  Dr. 
Wiesle,  the  A  matrix  was  rotated  from  the  nodal  frame  (nearly  an  inertial  frame)  into  the 
orbital  frame  (similar  to  the  CW  frame).  The  benefit  of  this  transfonnation  was  that 
many  of  the  elements  of  the  A  matrix  were  reduced  in  magnitude  resulting  in  a  maximum 
order  of  magnitude  variation  of  IE-3  over  the  period.  This  smaller  variation  arose 
speculation  that  the  time  varying  A  matrix  could  be  treated  as  a  constant  A  matrix. 

2.3.1  Nodal  Frame  A  Matrix.  The  original  dynamics  method  in  the  nodal  frame 
produced  an  A  matrix  with  element  variations  on  the  order  of  1E1  within  an  orbital 
period.  Explicit  element  variations  can  be  seen  in  Figure  6  below. 
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Nodal  Frame  A  Matrix  Elements 


Figure  6.  Variations  of  Each  A  Matrix  Element  Over  One  Orbit 


Most  notable  from  Figure  6  is  that,  while  a  majority  of  the  elements  are  explicitly 
equal  to  zero,  many  of  the  non-zero  elements  vary  periodically  with  different  amplitudes 
and  frequency  of  oscillation.  For  the  most  extreme  case,  the  maximum  amplitude 
difference  exceeded  one.  While  it  is  obvious  that  this  A  matrix  would  not  be  a  good 
candidate  for  a  linear  time  invariant  (LTI)  system,  this  assumption  was  used  for  the  nodal 
frame,  time-invariant  (AnTI)  model.  The  time-invariant  A  matrix  used  for  this  model 
was  taken  as  a  snap  shot  early  on  in  the  orbit  before  many  of  the  varying  elements  had 
the  opportunity  to  change 
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2.3.2  Orbital  Frame  A  Matrix.  In  an  effort  to  achieve  linear  time  invariance, 


the  original  dynamics  method  was  transformed  to  the  orbital  frame.  This  transformation 
provided  by  Dr.  Wiesel  produced  an  A  matrix  with  element  variations  on  the  order  of  1E- 
3  within  an  orbital  period.  Explicit  element  variations  can  be  seen  in  Figure  7  below. 


Orbital  Frame  A  Matrix  Elements 


Figure  7.  Variations  of  Each  A  Matrix  Element  Over  One  Orbit 


In  comparison  with  Figure  6,  it  is  clear  that  this  method,  while  still  producing 
some  periodically-varying  non-zero  elements  with  different  amplitudes  and  frequency  of 
oscillation,  the  amplitudes  are  much,  much  smaller  than  for  the  nodal  frame.  For  the 
most  extreme  case,  the  maximum  amplitude  difference  exceeded  IE-3  (see  Figure  9 
below).  This  A  matrix  appears  be  a  good  candidate  for  a  LTI  system,  so  this  assumption 
was  used  for  the  orbital  frame  A  matrix  time-invariant  (AoTI)  model.  The  time-invariant 
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A  matrix  used  for  this  model  was  also  taken  as  a  snap  shot  early  on  in  the  orbit  before  the 
varying  elements  had  the  opportunity  to  change 
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As  will  be  shown  in  the  results  chapters,  the  assumption  of  LTI  for  this  model 
worked  very  well  over  many  scenarios. 


2.3.3  Nodal  and  Orbital  Frame  Comparisons.  It  is  also  useful  to  look  at 
comparisons  of  different  elements  of  each  A  matrix. 

A  Matrix  Element  (4,1),  Max  Nodal  Variation 
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Figure  8.  A  Matrix  Element  (4,1)  Comparisons,  Max  Nodal  Variation 
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x  iq'3  A  Matrix  Element  (3,2),  Max  Orbital  Variation 


Figure  9.  A  Matrix  Element  (3,2)  Comparisons,  Max  Orbital  Variation 

x  iq’16  A  Matrix  Element  (1 ,1),  Typical  "Zero"  element 
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Figure  10.  A  Matrix  Element  (1,1)  Comparisons,  Typical  “Zero”  Element 


Many  of  the  zeros  for  the  orbital  frame  actually  varied  periodically  but  were  tiny. 
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III.  Estimation  and  Control  Methodology 


The  overriding  tasks  in  estimation  and  control  are  obtaining  a  good  estimate  for 
the  dynamic  states  and  balancing  the  desire  for  state  error  cancellation  (based  on  the  state 
estimates)  with  the  desire  to  keep  fuel  expenditure  very  low.  We  will  use  optimal  control 
methods,  specifically  Kalman  filter  and  linear  quadratic  regulator  theory,  to  strike  this 
balance.  In  order  to  complete  this  optimization  task,  we  must  first  set  up  the  conditions 
for  our  estimator  and  controller.  These  assumptions  and  rules  will  govern  the  building 
blocks  of  how  our  system  will  be  implemented. 

3.1  Finite  Time,  Non-linear  Time-varying  Problem 

Following  controller  design  theory  from  Bryson  [7],  we  start  with  the  more 
general  problem  where  the  state,  X  is  described  by 

X  =  f(X(t),U(t),t)  (24) 

This  handles  a  general  class  of  non-linear,  time-varying  systems.  A  well  known 
cost  function  for  these  problems  is  the  quadratic  perfonnance  index 

J  =  I .  Xj  •  Qf  •  Xf  +  i  •  |t‘f(XT  •  Q  ■  X  +  UT  •  R  ■  U)dt  (25) 

The  states  for  this  problem  are  typically  fonnulated  as  error  states  that  we  would 
like  to  regulate  as  close  to  zero  as  possible  for  a  given  amount  of  control  usage.  The 
general  form  of  the  problem  results  in  a  two-point  boundary  value  problem  which  cannot 
readily  be  solved  online.  To  make  our  lives  easier,  we  would  like  to  make  some 
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assumptions  to  give  us  a  solution  that  could  be  evaluated  online.  The  ideal  would  be  a 
solution  that  could  be  evaluated  on-board  a  satellite’s  navigation  computer. 

3.2  Finite  Time,  Linear  Time-invariant  Problem 

The  assumption  of  a  linear  time-invariant  system  allows  the  state  to  be  described 
by  the  state-space  equation: 

X  =  A-X  +  BU  (26) 

where  A  and  B  are  constant  matrices.  Using  a  quadratic  performance  mentioned 
in  the  previous  section,  the  optimal  control  solution  of  this  problem  is  given  by 

U  =  -K  •  X  (27) 

where  the  controller  gain  K  is  given  by 

K  =  R  -1  •  B  T  •  P(t)  (28) 

where  P(t)  is  the  solution  of  a  Riccati  differential  equation 

-  P  =  At  •  P(t)  +  P(t)  ■  A  -  P(t)  ■  B  ■  R_1  ■  Bt  ■  P(t)  +  Q  (29) 

While  this  Riccati  differential  equation  is  solvable,  the  process  would  not  be 
described  as  easy  or  useful  for  near  real-time  on-board  calculations.  Note  that  the  LTI 
system  has  an  optimal  control  law  with  time-varying  gains  [13]. 

3.3  Infinite  Horizon,  Linear  Time-invariant  Problem 

If  we  allow  our  final  time  to  approach  infinity,  we  can  assume  for  a  stable  system 
that  the  regulated  system  will  approach  steady  state.  As  tf  approaches  infinity,  the 
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performance  index  is  dominated  by  the  steady  state  nature,  so  P(t)  approaches  Pss.  As 
P(t)  approaches  a  constant  value,  P  approaches  zero.  Under  these  conditions,  the  Riccati 
differential  equation  from  the  previous  section  becomes  an  algebraic  Riccati  equation 

0  =  AT  •  Pss  +  Pss  •  A  -  Pss  •  B  •  R_1  •  BT  •  Pss  +  Q  (30) 

This  algebraic  Riccati  equation  is  readily  solvable,  and  this  relation  will  be  utilized 
in  Section  3.5.  However, we  have  made  many  assumptions  to  make  it  to  this  point,  and, 
while  time-invariance  would  be  nice  to  have,  the  states  for  our  orbiting  satellites  are  not 
time-invariant.  Due  to  the  repetitive  nature  of  orbiting  satellites  and  their  perturbations,  a 
time-varying  time  periodic  description  would  hold  more  true  (see  Section  2.1.2).  This 
time  periodic  assumption  will  be  used  in  the  next  section. 


3.4  Linear  Time  Periodic  Method  (LTP) 


The  importance  of  the  Riccati  equation  for  controls,  signals,  and  systems  has 
sparked  much  research  on  the  time-varying  periodic  Riccati  equation.  Following  a 
method  proposed  by  Bittanti,  et  al  [3],  the  following  state-space  equations  were  used  to 
describe  our  orbital  system: 


X(t)  =  A(t)  -X(t)  +  B  -U(t)  +  w(t) 
Y(t)  =  C  -X(t)  +  v(t) 


(31) 


Where  U(t)  was  the  control  input,  X(t)  was  the  state,  and  Y(t)  was  the  output. 
A(t)  was  a  periodic  matrix  of  period  T  such  that: 

A(t  +  T)  =  A(t)  (32) 

With  the  understanding  that  there  is  no  such  thing  as  a  completely  detenninistic 
process,  dynamics  noise  w  and  measurement  noise  v  were  used  to  describe  model 
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uncertainties,  non-linearities,  perturbations,  or  any  other  constraint  on  the  actual  system 
that  was  not  easily  described  mathematically.  As  will  be  tested  and  shown  later,  the 
measurement  noise  was  one  of  the  big  drivers  constraining  the  performance  of  the 
systems.  Variables  w  and  v  were  defined  as  uncorrelated  zero-mean  white  Gaussian 
noises  with  covariance 

E{w-wt}  =  VI  and  E{vvT}  =  V2  (33) 


B 


(34) 


For  this  thesis,  only  the  A  matrix  was  variable.  Both  the  B  and  C  matrices  were 
set  as  constant  matrices. 

fo  0  O^j 
0  0  0 
0  0  0 
1  0  0 
0  1  0 
V0  o  l) 

The  physical  meaning  of  this  B  matrix  was  that  for  all  time,  the  control  system 
had  perfect  firing  of  its  thrusters  in  each  direction  exclusively.  This  would  mean  that 
when  the  thruster  was  fired  in  the  y-direction,  all  the  thrust  was  concentrated  solely  in  the 
y-direction.  Further,  each  thruster  was  identical  in  perfonnance  capability.  This  was  an 
ideal  situation,  and  any  deviations  from  this  perfection  were  modeled  as  a  portion  of  the 
dynamics  noise  w  in  the  system.  The  standard  deviation  for  truth  model  dynamics  noise 
used  was  consistent  with  findings  from  a  paper  by  Bordner  and  Wiesel  that  a  comparable 
uncontrolled  satellite  would  drift  approximately  2  m/day  in  each  direction  [5]. 

r\  0  0  0  0  O^j 

C  =  0  1  0  0  0  0 

v0  0  1  0  0  0 ) 


(35) 
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The  physical  meaning  of  this  C  matrix  was  that,  in  the  absence  of  measurement 
noise  v,  the  control  system  would  receive  perfect  absolute  GPS  position  readings  in  the  x, 
y,  and  z  directions.  Again,  this  was  an  ideal  situation,  and  deviations  from  this  perfection 
were  modeled  as  the  measurement  noise  v  in  the  system.  Since  the  absolute  GPS 
position  error  doubles  in  the  vertical  direction,  the  standard  deviations  for  measurement 
noise  v  were  modeled  as  4  m  in  the  radial  x  direction  and  2  m  in  both  the  in  track  y 
direction  and  cross  track  z  directions  [4]. 

The  periodic  Riccati  equation  which  solves  the  state-space  system  of  equations  is 

-P(t)  =  A(t>P(t)+P(t>  A(t)T  -P(t)  C(t)T  V21  C(t)  P(t)+B(t)  Vl  B(t)T  (36) 


The  Hamiltonian  matrix  corresponding  to  the  periodic  Riccati  equation  is 


H(t)  = 


A(t)T 

-C(t)T  V2'1  C(t) 


-  B(t)  ■  VI  ■  B(t)T 
-  A(t) 


(37) 


This  Hamiltonian  matrix  has  symplectic  characteristics,  following  the  relationship 

H(t)  T  •  J  +  J-H(t)  =  0  (38) 

where  J  is  the  symplectic  matrix  consisting  of  zeros  and  identity  matrices 


J  = 


I 

0 


(39) 


The  Hamiltonian  transition  matrix  <J>H(t,t0)  also  has  symplectic  characteristics 
following  the  relationship, 

®H(t,t0)T-T<I>H(t,to)  =  J  (40) 

The  Hamiltonian  transition  matrix  can  be  partitioned  into  four  n  x  n  blocks, 


3>H  (t,to)  = 


<J>n(t,to)  d?12(t,to) 
d>2i(t,to)  d>22(t,to) 


(41) 


The  solution  to  the  periodic  Riccati  equation  with  initial  condition  Po  at  time  to  is 
P(t)  =  (<h21  (t,t0)  +  <h22(t,t0)  -Po)  ■  (  $ii  (t,t0)  +  <h12  (t,t0)  -Po)'1  (42) 
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Once  this  solution  was  computed  over  one  orbit,  it  could  be  used  to  compute  the 
optimal  control  gain  and  would  be  repeated  exactly  for  all  the  following  orbits.  Similar 
to  the  use  of  the  tenn  “exactly”  in  earlier  sections,  deviations  to  this  mode  were  modeled 
as  dynamics  noise  w  and  measurement  noise  v.  In  order  to  verify  that  this  solution  was 
accurate,  checks  of  the  symplectic  characteristics  of  each  of  the  building  blocks  were 
made.  For  details  on  these  checks,  their  anomalous  results,  and  the  many  attempts  to 
modify  the  system  to  get  positive  results,  see  Appendix  A. 

It  is  important  to  note  that  between  this  section  and  the  next,  there  is  an  important 
shift  in  both  the  underlying  method  and  the  tool  used.  As  this  section  attempted  to  solve 
the  periodic  Riccati  equation  using  code  written  in  the  fortran  computer  language,  the 
next  section  will  instead  solve  the  algebraic  Riccati  equation  using  the  standardized 
controls  routines  of  Matlab.  A  comparison  of  fortran  code’s  solution  of  the  algebraic 
Riccati  equation  to  Matlab ’s  solution  of  the  corresponding  algebraic  Riccati  equation 
yielded  errors  on  the  order  of  IE-29,  so  either  tool  could  handle  this  simpler  problem. 
Similar  to  Irvin’s  State  Dependent  Riccati  Equation  (SDRE)  method,  the  solution  to  the 
algebraic  Riccati  equation  was  assumed  valid  for  this  application  because,  in  effect,  we 
made  a  linear  solution  fit  a  curve  by  re-linearizing  for  each  of  the  25,000  time  steps  over 
the  orbit  [11]. 

3.5  Periodic  Riccati  Equation  -  Linear  Quadratic  Gaussian  (LQG)  Method 

Optimal  control  theory  provided  a  tunable  solution  to  our  periodic  problem  with 
the  linear  quadratic  Gaussian  (LQG)  method  combining  a  linear  quadratic  regulator  with 
a  Kalman  filter.  Following  Maciejowski,  the  state-space  form  of  our  problem  was  the 
same  as  in  the  previous  section  [17]. 
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X(t)  =  A(t)  -X(t)  +  B(t)  -U(t)  +  w(t) 


(43) 


Y(t)  =  C(t)  -X(t)  +  v(t) 

The  separation  principle  was  used  to  solve  this  LQG  problem  for  the  optimal 
solution.  Basically,  the  linear  quadratic  regulator  operates  on  X  ,  the  optimal  estimate  of 
state  X  from  the  Kalman  filter  much  as  it  would  if  the  full  state  were  available  as  a 
measurement  for  feedback.  By  combining  the  regulator  and  filter,  LQG  theory 
guarantees  at  least  stability  for  a  nominal  system,  but  not  necessarily  robustness. 

3.5.1  Kalman  Filter.  Kalman  filter  theory  provided  the  means  for  estimating  the 
state  X.  The  Kalman  filter,  as  seen  below  in  Figure  11,  would  take  inputs  of  the  plant 
input  U  and  plant  output  Y.  The  Kalman  filter  would  produce  an  output  of  the  state 
estimate  vector  X . 


Figure  1 1.  The  Kalman  Filter  [17] 


The  object  of  the  Kalman  filter  problem  is  to  minimize  the  cost 

J  =  E{[(X(t)  -  X(t)]  ■  W(t)  ■  [(X(t)  -  X(t)]T  }  (44) 

The  Kalman  filter  gain  matrix  Kf  was  calculated  as 

Kf=PfCT-N2'1  (45) 

where  Pf  was  the  solution  to  the  algebraic  Riccatti  equation 
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(46) 


Pf  At  +  A  Pf  -  Pf  CTN2'1-C  Pf  +  N1  =0 

The  filter  dynamics  noise  N 1  and  measurement  noise  N2  tuning  parameters  would 
be  varied  to  optimize  the  Kalman  filter’s  estimate  of  the  state  X  .  The  starting  point  was 
to  set  the  N1  value  to  the  expected  standard  deviation  squared  for  the  truth  model’s 
dynamics  noise  V 1  (standard  deviation  of  2  meters  per  day),  and  the  N2  value  was  varied 
as  orders  of  magnitude  of  the  expected  standard  deviation  squared  for  the  truth  model’s 
measurement  noise  V2  (standard  deviation  4  meters  in  the  radial  direction  and  2  meters 
in  the  in-track  and  cross-track  directions)  [4,  5].  Now  that  the  Kalman  filter  is  set  to  the 
“truth”  noises,  we  tune  the  filter  gain  by  varying  the  magnitude  of  N2  to  get  the  response 
we  want.  When  you  increase  the  magnitude  of  N2,  you  are  telling  the  filter  not  to  believe 
the  measurements  as  much.  When  you  decrease  the  magnitude  of  N2,  you  are  telling  the 
filter  that  the  measurements  it  receives  are  really  accurate,  and  it  should  believe  them 
more.  Usually,  it  is  better  to  tune  the  filter  first  for  one  noise  (N2),  and  then  the  other 
(Nl).  Measurement  noise  N2  was  chosen  as  the  first  tuner  because,  as  mentioned  earlier 
in  this  chapter,  measurement  noise  has  a  much  bigger  effect  on  system  performance  for 
our  system.  As  noted  later  in  the  thesis,  the  effects  of  changing  N 1  are  very  slight. 

3.5.2  Full  State  Feedback  Regulator.  The  object  of  the  regulator  problem  was 
to  come  up  with  the  control  signal  to  minimize  the  cost 

T 

J  =  J ((Xdes  -  X)T  •  Q  ■  (Xdes  -  X)  +  UT  •  R  ■  U)dt  (47) 

0 

from  the  initial  state-space  description  in  Equation  43.  The  solution  was  to  keep 
control  signal  U  as  a  linear  function  of  the  state,  or  in  this  case  of  the  state  error  (Xdes  -  X) 

U  =  -  Kc-  (Xdes  -  X)  (48) 

The  optimal  state-feedback  matrix  Kc  was  calculated  as 
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(49) 


Kc  =  R'1-Bt-Pc 

where  Pc  was  the  solution  to  the  algebraic  Riccati  equation 

At-Pc  +  Pc- A  -  Pc-B-R'1  -Bt-Pc  +  Q  =  0  (50) 


Figure  12.  Optimal  Full  State  Feedback  Regulator 


Attached  to  the  truth  model,  the  optimal  full  state  feedback  regulator  would  try  to 
drive  the  state  error  (Xdes  -  X)  to  zero.  The  Q  and  R  values  were  tuned  to  produce  the 
optimal  response  via  gain  Kc.  Following  a  method  proposed  by  Irvin,  the  Q  value  was 
initially  set  to  identity,  and  then  the  R  values  were  varied  by  orders  of  magnitude  [9]. 
Once  the  acceptable  R  value,  in  tenns  of  control  usage  and  RMS  error,  was  found,  Q 
values  were  tweaked  to  see  if  improvements  could  be  made.  The  process  is  explained 
more  in  the  following  chapters,  but  eventually  controller  tuning  values  of  R  =  IE-4  and  Q 
=  1E0  were  found  to  be  the  best  compromise  for  comparisons  of  the  different  models. 

The  truth  model  was  used  to  output  not  only  the  measurements  of  the  true  state  as 
calculated  by  the  fortran  routines  from  Section  2.2.1,  but  it  also  added  the  perturbing 
noises  w  and  v  (of  variance  V 1  and  V2)  to  account  for  un-modeled  disturbances.  Note 
that  the  truth  model  was  not  changed  during  the  tuning.  In  essence,  we  cannot  change 
truth,  but  by  changing  the  gain  of  the  regulator,  we  controlled  how  the  regulator 
interpreted  what  it  saw  (Xdes  -  X),  and,  thus,  how  it  reacted  (U). 
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3.5.3  Output  Feedback  Compensator.  The  Kalman  filter  and  controller  were 
married  together  in  series  to  form  a  new  optimal  output  feedback  compensator  as  seen  in 
Figure  13.  Similar  to  the  regulator  tuning,  variations  on  Q  and  R  were  used  to  optimize 
the  controller  gain  Kc.  To  provide  even  better  results,  the  filter  N1  and  N2  tuning  values 
were  tweaked  as  in  Section  3.5.1  to  optimize  the  Kalman  filter  gain  Kf.  Additional 
details  on  the  actual  tuning  appears  in  Appendix  B. 

Once  the  model  was  tuned,  additional  scenarios  involving  the  different  A  matrices 
from  Section  2.2.2  were  executed  to  compare  the  degradation  when  switching  from  time 
varying  control  to  time  invariant  control.  Note  that  when  we  switch  from  time-varying  to 
time-invariant  control,  the  truth  model  must  remain  the  same. 


Figure  13.  Output  Feedback  Compensator  [17] 


32 


IV.  Full  State  Feedback  Regulator  Results 


In  order  for  the  regulator  and  compensator  to  be  tuned  and  checked  for  validity, 
Monte  Carlo  data  runs  were  executed.  For  response  to  initial  conditions,  only  one  Monte 
Carlo  run  was  accomplished;  multiple  runs  would  be  needed  for  statistical  measurements. 
For  both  the  regulator  and  compensator,  the  scenarios  were  originally  run  with  an  initial 
perturbation  of  10  meters  in  the  radial  direction.  Some  later  runs  used  initial 
perturbations  of  greater  magnitudes  or  in  other  directions,  but  for  comparative  purposes, 
the  10  meter  initial  radial  perturbation  was  used  as  the  standard.  This  perturbation  would 
be  consistent  with  the  perturbation  caused  by  increased  solar  activity  such  as  the 
increased  air  drag  due  to  a  small  solar  flare.  This  chapter  will  present  the  results  with 
figures  and  tables  focused  on  the  control  and  RMS  error.  The  RMS  error  was  calculated 
as  steady  state  RMS  error  by  averaging  the  RMS  error  over  one  orbit  with  no  initial 
perturbation.  Settling  times  were  figured  using  a  5%  of  steady  state  threshold. 

4.1  Full  State  Feedback  Regulator,  Initial  Perturbation  Dx  =  10  meters 

The  first  step  of  the  tuning  process  for  the  output  feedback  compensator  was  to 
set  up  and  tune  the  full  state  feedback  regulator.  As  mentioned  in  the  previous  chapter, 
the  full  state  feedback  regulator  is  more  of  an  idealized  system  where  your  controller  can 
see  all  of  the  states  without  any  measurement  noise  v  distorting  the  measured  feedback. 
As  can  be  seen  in  Figure  14  below,  although  the  system  does  not  receive  measurement 
noise  v  in  the  feedback,  the  dynamics  noise  w  of  the  truth  model  does  factor  into  the 
system.  But  as  noted  in  the  last  chapter,  the  magnitudes  of  the  dynamics  noise  (standard 
deviation  of  2  m/day)  were  small  enough  for  the  system  to  handle  them  without  much 
problem. 
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Although  ideal,  the  results  of  the  regulator  tuning  process  proved  to  be  useful  in 
the  later  tuning  of  the  more  realistic  compensator.  Results  of  control  versus  settling  time 
and  control  versus  RMS  error  were  compared  to  find  the  optimum  values  for  the  R  tuning 
parameter.  Although  data  will  be  presented  for  all  the  tuning  values,  comments  will 
focus  on  results  for  the  tuning  value  R=lE-4. 


4.1.1  Time-invariant  Anodal  Model  (AnTI).  The  first  model  tested  was  the 
case  of  the  time-invariant  A  matrix  in  the  nodal  frame.  The  physical  meaning  of  this 
model  is  that  the  controller  receives  full  state  information  and  controls  the  system  based 
on  a  time-invariant  model  (single  representative  snapshot)  of  the  dynamics  of  the  system. 


Figure  14.  Full  State  Feedback  Regulator,  AnTI 


Due  to  the  variances  in  the  components  of  the  nodal  frame  A  matrix  as 
detailed  in  Section  2.3,  this  method  is  not  expected  to  perfonn  as  well  as  the  other 
methods.  Surprisingly,  the  idealistic  nature  of  the  regulator  provides  adequate  control  for 
this  scenario. 
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Regulator  Error  (x  -  xdes),  Dx=10m,  R=1E-4,  AnTI 


Figure  15.  Full  State  Feedback  Regulator  Error  for  R=lE-4,  AnTI 


Regulator  Control,  Dx=10m,  R=1E-4  ,  AnTI 


Figure  16.  Full  State  Feedback  Regulator  Control  for  R=lE-4,  AnTI 
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4.1.2  Time-varying  Anodal  Model  (AnTV).  The  time-varying  A  matrix  in  the 
nodal  frame  was  the  next  regulator  model  tested.  For  this  model  the  controller  receives 
the  full  state  information  and  controls  the  system  knowing  the  time  periodic  dynamics  of 
the  system  at  each  step  of  the  orbit. 


Truth  Model 
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Figure  17.  Full  State  Feedback  Regulator,  AnTV 


The  representative  plots  for  error  and  control  over  one  orbit  show  acceptable 
response  to  the  10  meter  initial  perturbation  for  the  selected  tuning  value  R=lE-4. 

In  comparison  of  Figure  16  and  Figure  19,  the  control  responses  for  the  time- 
varying  and  time-invariant  control  are  very  similar.  This  shows  that  for  this  ideal  set-up 
where  the  controller  was  able  to  see  the  exact  states  without  measurement  noise,  even  the 
models  that  are  not  well  described  by  LTI,  can  produce  acceptable  responses. 
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Regulator  Error  (x  -  xdes),  Dx=10m,  R=1E-4,  AnTI 


Figure  18.  Full  State  Feedback  Regulator  Error  for  R=lE-4,  AnTV 


Regulator  Control,  Dx=10m,  R=1E-4,  AnTV  (detail) 


Figure  19.  Full  State  Feedback  Regulator  Control  for  R=lE-4,  AnTV 
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4.1.3  Time-varying  Aorbital  Model  (AoTV).  The  third  model  tested  was  the 
case  of  the  time -varying  A  matrix  in  the  orbital  frame.  Similar  to  AnTV,  the  physical 
meaning  of  this  model  is  that  the  controller  receives  the  full  state  information  and 
controls  the  system  knowing  the  dynamics  of  the  system  at  each  step  of  the  orbit.  This 
time  the  dynamics  of  the  system  are  in  the  orbital  frame.  Because  of  the  new  orbital 
frame  of  reference  and  the  inherent  adaptability  of  the  time-varying  model,  this  method  is 
expected  to  provide  the  best  performance. 


Figure  20.  Full  State  Feedback  Regulator,  AoTV 


The  error  plot  in  Figure  21  shows  how,  due  to  orbital  dynamics,  a  perturbation  in 
the  radial  direction  will  cause  the  satellite  to  slow  down,  and  thus  degrade  the  y  in-track 
direction.  This  type  of  cross-correlation  would  need  further  study,  and  adjustments  to  the 
controller  gain  would  be  an  area  for  improvement  for  the  future. 
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Regulator  Error  (x  -  xdes),  Dx=10m,  R=1E-4,  AoTV 


Figure  2 1 .  Full  State  Feedback  Regulator  Error  for  R=lE-4,  AoTV 


Regulator  Control,  Dx=10m,  R=1E-4,  AoTV  (detail) 


Figure  22.  Full  State  Feedback  Regulator  Control  for  R=lE-4,  AoTV 
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4.1.4  Time-invariant  Aorbital  Model  (AoTI).  The  time-invariant  A  matrix  in 
the  orbital  frame  was  the  last  regulator  model  tested.  For  this  model,  the  controller 
receives  the  full  state  information  and  controls  the  system  knowing  only  a  single  time- 
invariant  snapshot  of  the  dynamics  of  the  system  at  each  step  of  the  orbit.  This  time  the 
dynamics  of  the  system  are  in  the  orbital  frame. 


Figure  23.  Full  State  Feedback  Regulator,  AoTI 


Regulator  Error  (x  -  xdes),  Dx=10m,  R=1E-4,  AoTI 


Figure  24.  Full  State  Feedback  Regulator  Error  for  R=lE-4,  AoTI 
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Figure  25.  Full  State  Feedback  Regulator  Control  for  R=lE-4,  AoTI 


4.1.5  Regulator  Model  Comparisons.  The  results  for  the  four  models  were 
compared  for  control,  RMS  error,  and  settling  time  over  an  order  of  magnitude  range  of 
controller  tuning  values  from  R=lE-6  to  R=1E0.  Each  data  point  on  the  figures  below 
represents  the  value  for  one  particular  tuning  value  for  R.  For  completeness,  R  values 
above  and  below  this  range  were  tested  to  ensure  the  trends  of  the  data  did  not  change  for 
the  extremes. 
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Control  (meters/second) 


Total  Regulator  Control  Over  One  Orbit,  Dx=10m 


Figure  26.  Full  State  Feedback  Regulator  Control,  Tuning  R’s 


Regulator  Settling  Times  Over  One  Orbit,  Dx=10m 


Figure  27.  Full  State  Feedback  Regulator  Settling  Times,  Tuning  R’s 
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For  three  out  of  four  models,  the  results  of  the  control  versus  settling  time  plot 
clearly  favored  a  tuning  value  of  R=lE-4.  The  one  undecided  model,  AnTI,  was  not  sure 
if  it  liked  R=lE-4  (which  gave  it  lower  control)  or  R=lE-3  or  R=lE-2  (which  gave  it 
lower  settling  time).  Physically,  the  R=lE-4  value  translates  to  telling  the  controller  that 
the  errors  it  should  correct  are  reliable  but  not  to  be  trusted  fully.  The  goal  of  the  tuning 
process  is  to  minimize  both  control  and  settling  time,  which  corresponds  to  the  lower  left 
comer  of  Figure  28  below.  Note  that  the  lower  magnitudes  of  tuning  value  R  correspond 
to  the  leftmost  data  points  on  the  figure  below. 


Regulator  Settling  Time  vs  Control  Over  One  Orbit,  Dx=10m  (detail) 


Figure  28.  Regulator  Settling  Time  vs.  Control,  Tuning  R’s 


An  important  result  to  note  is  the  consistency  of  the  time-varying  and  time- 
invariant  results  for  the  orbital  A  matrix.  As  can  be  seen  in  the  figure  above  as  well  as 
Tables  2-4,  control,  RMS  errors,  and  settling  times  were  nearly  identical  for  AoTI  and 
AoTV  for  all  values  for  R.  Referring  back  to  Section  2.3.2,  where  the  A  matrix 
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components  in  the  orbital  frame  are  compared  over  one  orbit,  it  becomes  clear  why  AoTI 
and  AoTV  should  produce  very  similar  results.  That  was  the  reason  for  this  frame 
change:  to  get  an  A  matrix  that  looks  very  close  to  being  constant  over  the  whole  orbit. 
The  variability  of  the  nodal  frame  components  (see  Section  2.3.1)  could  be  a  reason  as  to 
why  the  orbital  frame  models  consistently  required  14%  less  fuel.  The  complexity  of  the 
model  trying  to  manipulate  all  the  widely  varying  periodic  changes  of  the  dynamics  could 
cause  the  system  to  overwork  itself,  spending  more  fuel.  The  nodal  controller  might 
overcompensate  and  under-compensate  more  compared  to  the  orbital  frame  models, 
where  the  system  knows  pretty  well  what  its  baseline  is.  Notice  the  big  overshoots  for 
the  control  plots  for  the  nodal  frame  models  versus  the  small  overshoots  for  the  orbital 
frame  models.  Although  the  time-invariant  Anodal  scenario  varied  a  little  from  the  rest, 
for  the  target  tuning  values  of  IE-6,  IE-5,  and  IE-4,  time-varying  and  time-invariant 
methods  produced  control,  RMS  errors,  and  settling  times  nearly  the  same  for  AnTV  and 
AnTI. 


able  2.  Regulator  Total  Control  Over  One  Orbit,  Dx  =  10  m 


Control  (m/sec) 

R=lE-6 

R=lE-5 

R=lE-4 

R=lE-3 

R=lE-2 

R=1E-1 

R=1E0 

AnTI 

0.0343 

0.0331 

0.0329 

0.0335 

0.0359 

0.0459 

0.0969 

AnTV 

0.0343 

0.0330 

0.0329 

0.0335 

0.0356 

0.0394 

0.0449 

AoTI 

0.0304 

0.0291 

0.0289 

0.0293 

0.0304 

0.0315 

0.0323 

AoTV 

0.0304 

0.0291 

0.0289 

0.0293 

0.0304 

0.0315 

0.0323 
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Table  3.  Regulator  Settling  Times,  Dx  =  10  m 


(% 

R=lE-6 

R=lE-5 

R=lE-4 

R=lE-3 

R=lE-2 

R=1E-1 

R=1E0 

orbit) 

AnTI 

41.31 

41.22 

40.94 

40.04 

37.48 

50.31 

>100 

AnTV 

41.36 

41.39 

41.48 

41.75 

42.39 

41.86 

34.56 

AoTI 

41.37 

41.40 

41.49 

41.77 

42.48 

43.52 

52.86 

AoTV 

41.37 

41.40 

41.50 

41.77 

42.48 

43.53 

52.86 

Table  4.  Regulator  Average  RMS  Error  Over  One  Orbit,  Dx  =  0m 


Error 

R=lE-6 

R=lE-5 

R=lE-4 

R=lE-3 

R=lE-2 

R=1E-1 

R=1E0 

(meters) 

AnTI 

0.79E-4 

0.54E-4 

0.57E-4 

0.78E-4 

0.52E-4 

0.65E-4 

0.62E-4 

AnTV 

0.49E-4 

0.67E-4 

0.75E-4 

0.48E-4 

0.86E-4 

0.93E-4 

0.61E-4 

AoTI 

0.59E-4 

0.60E-4 

0.96E-4 

0.51E-4 

0.60E-4 

0.77E-4 

0.51E-4 

AoTV 

0.77E-4 

0.68E-4 

0.70E-4 

0.56E-4 

0.45E-4 

0.87E-4 

0.87E-4 
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V.  Output  Feedback  Compensator  Results 


After  tuning  the  full  state  feedback  regulator,  the  next  step  was  to  tune  the  output 
feedback  compensator.  As  mentioned  in  previous  chapters,  the  output  feedback 
compensator  is  a  more  realistic  system  where  the  compensator  can  only  see  the  noisy 
measurement  of  the  position  states.  Comparing  the  compensator  control  plots  to  the 
corresponding  regulator  ones,  it  is  easily  seen  how  the  noise  degrades  the  system.  But 
the  results  of  the  more  realistic  compensator  proved  similar  to  those  of  the  idealized 
regulator.  The  AnTI  model  was  the  exception.  It  diverged  for  all  scenarios,  and  thus,  it 
will  be  shown  in  Section  5.1,  but  will  not  be  presented  in  the  later  comparison  scenarios. 
As  mentioned  in  the  previous  chapter,  RMS  error  was  taken  as  a  steady  state  value.  RMS 
error  tested  each  specific  noise  level  or  step  size  condition  to  the  same  truth  model  with 
no  initial  perturbations.  Results  of  control  versus  settling  time  and  RMS  error  were 
compared  to  find  the  optimum  values  for  the  controller  Q  and  R  tuning  parameters  as 
well  as  the  filter  N1  and  N2  tuning  parameters. 

5.1  Output  Feedback  Compensator,  Initial  Perturbation  Dx  =  10  meters 

As  mentioned  in  the  previous  section,  this  scenario,  with  an  initial  perturbation  of 
10  meters  in  the  radial  direction,  was  the  standard  that  will  be  used  to  compare  to  for  all 
the  scenarios  in  following  sections. 


5.1.1  Time-invariant  Anodal  Model  (AnTI). 

case  of  the  time-invariant  A  matrix  in  the  nodal  frame. 


The  first  model  tested  was  the 
This  is  similar  to  the  AnTI 
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regulator  model.  But  this  time  the  controller  receives  the  noisy  output  information  and 
controls  the  system  with  only  a  time-invariant  snapshot  of  the  dynamics  of  the  system. 


Compensator  Error  (x  -  xdes),  Dx=10m,  R=1  E-4,  AnTI 


Figure  30.  Output  Feedback  Compensator  Error  for  R=lE-4,  AnTI 
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Figure  3 1 .  Output  Feedback  Compensator  Control  for  R=lE-4,  AnTI 


The  AnTI  position  errors  for  all  tuning  values  for  R  would  be  classified  as  slowly 
diverging  over  the  single  orbit.  Continuing  past  the  one  first  orbit  (see  Figure  32  below), 
the  AnTI  model  continues  to  diverge  with  no  hope  of  recovery  for  any  tuning  value  for  R. 
This  trend  held  true  for  all  of  the  compensator  scenarios  tested. 
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Compensator  Error  (x  -  xdes),  Dx=10m,  R=1E-4,  AnTI 


Figure  32.  Compensator  Error  Over  Three  Orbits,  Dx=10m,  AnTI 


5.1.2  Time-varying  Anodal  Model  (AnTV).  The  next  compensator  model 
tested  was  the  case  using  the  time-varying  A  matrix  in  the  nodal  frame.  The  physical 
meaning  of  this  model  is  that  the  controller  receives  the  noisy  output  information  and 
controls  the  system  knowing  the  dynamics  of  the  system  at  each  step  of  the  orbit. 
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Control  (meters/second) 


Compensator  Error  (x  -  xdes),  Dx=10m,  R=1  E-4,  AnTV 


Figure  34.  Output  Feedback  Compensator  Error  for  R=lE-4,  AnTV 


Compensator  Control,  Dx=10m,  R=1  E-4,  AnTV  (detail) 
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Figure  35.  Output  Feedback  Compensator  Control  for  R=lE-4,  AnTV 
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5.1.3  Time-varying  Aorbital  Model  (AoTV).  The  time-varying  A  matrix  in  the 
orbital  frame  was  the  third  model  tested.  For  this  model,  the  controller  receives  the  noisy 
output  information  and  controls  the  system  knowing  the  dynamics  of  the  system  at  each 
step  of  the  orbit.  This  time  the  dynamics  of  the  system  are  in  the  orbital  frame. 


Compensator  Error  (x  -  xdes),  Dx=10m,  R=1  E-4,  AoTV 


Figure  37.  Output  Feedback  Compensator  Error  for  R=lE-4,  AoTV 
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Figure  38.  Output  Feedback  Compensator  Control  for  R=lE-4,  AoTV 


5.1.4  Time-invariant  Aorbital  Model  (AoTI).  The  final  compensator  model 
tested  was  the  time-invariant  orbital  frame  A  matrix.  Similar  to  AnTI,  the  physical 
meaning  of  this  model  is  that  the  filter  and  the  controller  both  receive  the  noisy  output 
information  and  control  the  system  with  only  a  single  representative  snapshot  of  the 
dynamics  of  the  system.  This  time  the  dynamics  of  the  system  are  in  the  orbital  frame. 
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Compensator  Error  (x  -  xdes),  Dx=10m,  R=1E-4,  AoTI 


Figure  40.  Output  Feedback  Compensator  Error  for  R=lE-4,  AoTI 


Figure  41 .  Output  Feedback  Compensator  Control  for  R=lE-4,  AoTI 
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5.1.5  Compensator  Model  Comparisons.  The  results  for  the  three  models  were 
compared  for  control,  RMS  error,  and  settling  time  over  an  order  of  magnitude  range  of 
controller  tuning  values  from  R=lE-6  to  R=1E0.  This  range  of  magnitudes  was  validated 
through  the  full-state  feedback  regulator  tuning.  To  ensure  completeness,  tuning  R 
values  above  and  below  this  range  were  also  tested,  but  not  shown,  to  ensure  the  trends  of 
the  data  did  not  change  for  the  extremes  of  the  output  feedback  compensator  as  well. 
Figures  for  the  converging  models  AnTV,  AoTI,  and  AoTV  are  presented  first  with 
comparative  plots  for  the  diverging  model  AnTI  following  after. 

Note  that  in  the  previous  chapter  for  the  regulator  control  plots,  we  saw  the  nodal 
frame  models  with  bigger  overshoots  than  the  orbital  frame  models.  Now  that 
measurement  noise  is  added,  we  see  the  same  basic  control  plots  (a  little  noisier)  for  the 
compensator  orbital  frame  models,  but  now  for  the  compensator  the  nodal  frame  models 
really  start  to  show  how  the  noise  corrupts  their  attempts  to  control  the  system.  Compare 
Figure  6  from  Section  2.3.1  with  Figure  35  from  Section  5.1.2.  It  becomes  noticeable 
that,  as  the  components  of  the  A  matrix  begin  to  periodically  propagate  from  their  initial 
values,  the  noise  drives  the  control  to  become  more  variable.  It  is  important  to  note  here 
what  the  A  components  actually  represent  are  changes  to  the  dynamics  state  over  each 
step.  This  A  matrix  is  not  diagonal,  so  the  component  magnitudes  represent  correlations 
between  positions  and  velocities  that  increase  and  decrease  throughout  the  orbit.  From 
Equation  43,  the  total  change  in  the  states  includes  the  dynamics  state  terms  plus  control 
plus  dynamics  noise.  Since  the  expected  value  of  the  noise  at  different  points  of  the  orbit 
stays  the  same,  the  noise  will  play  a  bigger  role  when  the  A  matrix  components  are 
smaller  such  as  at  the  initial  time.  But  as  time  approaches  20-40%  of  the  orbit,  control 
quiets  down  because  the  A  matrix  components  periodically  increase  their  magnitudes  as 
they  propagate  through  their  orbits.  The  A  matrix  dynamics  start  to  dominate  over  the 
noise,  thus  you  see  less  variation  in  the  control  plots.  But  once  again,  as  the  A  matrix 
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components  start  to  decrease  towards  their  initial  values,  the  noise  again  starts  to 
dominate  and  drives  the  control  to  be  more  variable.  Looking  at  the  A  matrix  elements, 
most  of  them  propagate  to  zero  at  the  1/2  orbit  point  of  time.  It  is  not  as  noticeable  in 
these  plots,  but  due  to  its  periodic  nature,  the  orbital  frame  models  should  also  follow  this 
trend  to  a  much  lesser  degree. 


Total  Compensator  Control  Over  One  Orbit,  Dx=10m 


Figure  42.  Compensator  Control,  Dx=10m,  Tuning  R’s 
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Compensator  Average  RMS  Errors  Over  One  Orbit,  Dx=0m  (over  10  runs) 


Figure  43.  Compensator  RMS  Error,  Dx=10m,  Tuning  R’s 


Compensator  Settling  Times  Over  One  Orbit,  Dx=10m 


Figure  44.  Compensator  Settling  Times,  Dx=10m,  Tuning  R’s 
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Similar  to  the  regulator  results,  the  above  RMS  error  data  (Figure  43)  is  not  very 
helpful  for  using  trends  to  tune  specific  R  values,  but  it  does  show  a  trend  between  the 
orbital  frame  models  and  the  nodal  frame  model.  The  nodal  frame  model  consistently 
exceeds  an  order  of  magnitude  over  the  orbital  frame  models.  This  trend  ties  into  the 
dynamics  and  noise  correlation  discussed  earlier  in  this  section. 

For  settling  times,  Figure  44,  the  AnTV  value  drops  dramatically  for  the  R=1E0 
tuning  value.  This  trend  is  the  same  as  was  seen  with  the  regulator  in  Figure  27.  Both  of 
the  nodal  frame  regulator  models  had  a  minimum  settling  time  value  one  or  two  tuning  R 
magnitudes  below  the  tuning  value  for  which  they  would  diverge.  For  the  regulator, 
AnTI  minimized  settling  time  at  R=lE-2,  while  AnTV  was  minimized  at  R=1E0  similar 
to  the  compensator.  For  the  AnTV  compensator  model,  the  next  tuning  R  value  of 
R=1E1  drove  the  model  to  diverge.  This  was  a  feature  of  the  nodal  frame  models,  and 
not  for  the  orbital  frame  models.  The  orbital  frame  models  were  tested  over  many 
intennediate  R  values,  but  they  produced  no  minimizing  settling  time  dip  before  they 
reached  the  R  value  where  the  model  diverged.  The  reason  behind  the  minimum  settling 
time  is  tied  to  overshoot.  The  AnTV  model  in  Figure  34  is  over-damped  with  no 
overshoot,  so  it  takes  it  a  while  to  reach  the  5%  steady  state  value  that  settling  time 
requires.  The  minimized  settling  times  occur  when  the  system  is  under-damped  when  the 
system  quickly  overshoots  and  then  settles  into  the  5%  steady  state  value  threshold.  The 
next  magnitude  for  tuning  R  usually  provides  a  diverging  system  because  the  overshoot 
becomes  too  large  for  the  system  to  recover  from.  The  orbital  frame  models  did  not 
produce  similar  trends  because,  for  all  presented  R  values,  the  error  was  slightly  under¬ 
damped.  For  larger  R  values,  the  orbital  frame  models  gracefully  became  more  under¬ 
damped.  In  most  cases,  an  R  value  of  1E1  drove  most  systems  to  divergence.  Since  this 
feature  only  minimized  settling  time  for  nodal  frame  models  and  did  not  reduce  RMS 
error  or  control,  it  was  noted,  but  not  investigated  further. 
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Compensator  Settling  Time  vs  Control  Over  One  Orbit,  Dx=10m 


Settle  Time  (Orbits) 

Figure  45.  Compensator  Settling  Time  vs.  Control,  Dx=10m 


As  before  with  the  regulator,  the  best  balance  between  control  and  settling  time 
for  the  compensator  was  R=lE-4.  For  both  orbital  frame  models,  the  tuning  R  value  for 
best  control  was  R=lE-4,  but  AnTV  minimized  control  for  R=lE-2.  While  all  the 
models  reduced  settling  time  for  smaller  R  values,  for  comparative  purposes  a  single 
value  R=lE-4  provided  an  acceptable  compromise  between  minimal  control  and  minimal 
settling  time.  Although  the  consistency  of  the  time-varying  and  time-invariant  results  for 
the  orbital  frame  A  matrix  are  not  the  same  as  with  the  regulator,  they  are  still  pretty 
close.  It  is  also  a  pleasant  surprise  to  see  that  the  compensator’s  control  and  settling  time 
results  for  AoTI  and  AoTV  deviated  only  slightly  over  the  regulator  results.  Comparing 
regulator  to  compensator  results  for  both  AoTI  and  AoTV  (see  Table  5  below),  control 
degraded  less  than  one  percent.  Even  with  the  handicap  of  noisy  position  data,  the 
compensator,  with  its  combined  filter  and  controller,  was  able  to  control  the  system 
almost  as  well  as  an  idealistic  full  state  feedback  regulator,  but  the  RMS  error  suffered. 
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Table  5.  Compensator/Regulator  Comparisons  for  Dx=10m,  R=TE-4 


Total  Control  (1  orbit) 

Settling  Time 

RMS  error  (Dx=0m) 

(m/sec) 

%  over  Reg 

(%  orbit) 

%  over  Reg 

(meters) 

%  over  Reg 

AnTI 

0.3431  * 

90.4  %  * 

>  100  * 

59.1  % 

0.1770 

100.0  % 

AnTV 

0.0471 

30.3  % 

48.1 

-13.8 

0.0175 

99.6  % 

AoTI 

0.0291 

0.7  % 

41.0 

1.1  % 

0.0012 

91.9% 

AoTV 

0.0291 

0.8  % 

41.4 

0.2  % 

0.0015 

95.3  % 

*  Diverging  System 


Total  Compensator  Control  Over  One  Orbit,  Dx=10m 


Figure  46.  Compensator  Control,  Dx=10m,  Tuning  R’s 


The  AnTI  model  followed  the  same  trend  as  the  others  over  the  range  of  R  values, 
but  even  the  best  tuning  for  control,  R=lE-4,  yielded  control  results  which  were  an  order 
of  magnitude  greater  than  those  for  the  other  scenarios.  The  order  of  magnitude 
difference  is  not  surprising  considering  the  diverging  nature  of  the  results.  The  trend  is 
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important,  because  it  shows  that  the  model  is  trying  to  control  the  system,  but  the  time- 
invariant  model  we  gave  it  was  inadequate.  Since  this  model  diverged  for  all  tuning  R 
values  and  for  all  compensator  scenarios,  its  results  are  not  included  in  the  later  sections. 


Table  6.  Compensator  Total  Control  Over  One  Orbit,  Dx  =  10  m,  Tuning  R’s 


Control  (m/sec) 

R=lE-6 

R=lE-5 

R=lE-4 

R=lE-3 

R=lE-2 

R=1E-1 

R=1E0 

AnTI 

0.354  * 

0.344  * 

0.343  * 

0.351  * 

0.380  * 

0.442  * 

0.486  * 

AnTV 

0.0907 

0.0617 

0.0471 

0.0406 

0.0388 

0.0410 

0.0461 

AoTI 

0.0315 

0.0296 

0.0291 

0.0294 

0.0304 

0.0315 

0.0323 

AoTV 

0.0314 

0.0296 

0.0291 

0.0295 

0.0304 

0.0315 

0.0323 

*  Diverging  system 


Table  7.  Regulator  Settling  Times,  Dx  = 


0  m,  Tuning  R’s 


(%  orbit) 

R=lE-6 

R=lE-5 

R=lE-4 

R=lE-3 

R=lE-2 

R=1E-1 

R=1E0 

AnTI 

>  100  * 

>  100  * 

>  100  * 

>  100  * 

>  100  * 

>  100  * 

>100  * 

AnTV 

47.98 

48.01 

48.11 

48.41 

49.20 

49.08 

37.08 

AoTI 

40.92 

40.95 

41.03 

41.29 

41.98 

43.06 

52.56 

AoTV 

41.28 

41.31 

41.40 

41.68 

42.40 

43.47 

53.16 

*  Diverging  system 


Table  8.  Compensator  Average  RMS  Error  Over  One  Orbit,  Dx  =  0  m,  Tuning  R’s 


Error  (meters) 

R=lE-6 

R=lE-5 

R=lE-4 

R=lE-3 

R=lE-2 

R=1E-1 

R=1E0 

AnTI 

0.179 

0.179 

0.177 

0.171 

0.148 

0.071 

0.098 

AnTV 

0.0367 

0.0120 

0.0175 

0.0445 

0.0398 

0.0218 

0.0354 

AoTI 

0.0012 

0.0014 

0.0012 

0.0035 

0.0014 

0.0033 

0.0024 

AoTV 

0.0013 

0.0016 

0.0015 

0.0016 

0.0009 

0.0017 

0.0010 
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5.1.6  Filter  Tuning  N1  and  N2.  With  the  controller  properly  tuned  at  R=lE-4, 
the  next  step  was  to  tweak  the  filter  tuning  values  N1  and  N2.  Following  a  similar 
method  as  with  the  Q  and  R,  the  N 1  value  was  set  to  the  expected  standard  deviation 
squared  for  the  truth  model’s  dynamics  noise  (standard  deviation  of  2  meters  per  day), 
and  the  N2  value  was  varied  as  orders  of  magnitude  of  the  expected  standard  deviation 
squared  for  the  truth  model’s  measurement  noise  (standard  deviation  of  4  meters  in  the 
radial  direction  and  2  meters  in  the  in- track  and  cross-track  directions).  Again,  there  was 
a  trade-off  between  total  control  used  and  average  steady  state  RMS  error.  Tables  and 
figures  of  comparative  results  for  various  filter  N2  tuning  values  for  the  favored 
controller  tuning  value  R=lE-4  are  listed  in  Appendix  B.  The  most  notable  result  of  this 
tuning  process  was  that  either  the  control  or  RMS  error  could  be  improved  over  the 
nominal  N2=1E0,  but  not  both.  Steady  state  RMS  error  is  reduced  by  increasing  N2,  in 
effect,  telling  the  filter  to  believe  its  dynamics  model  better  than  the  noisy  measurements. 
Once  you  have  a  perturbation  to  deal  with,  the  control  and  settling  time  favor  the  smaller 
N2  values,  placing  more  emphasis  on  the  noisy  measurements  over  the  dynamics  model. 
Each  model  preferred  a  different  magnitude  of  N2.  Although  a  slight  amount  of  benefit 
could  be  gained  by  more  rigorously  tuning  each  filter,  for  comparative  purposes,  the  filter 
tuning  value  was  set  to  the  truth  model  variances  without  increasing  or  decreasing  the 
tuning  value. 


5.2  Initial  Perturbation  Dy  =  10  meters 

The  previous  section’s  initial  perturbation  of  10  meters  was  in  the  radial  direction, 
which  should  be  the  worst  direction  for  the  models  since  the  accuracy  of  GPS  readings  is 
typically  halved  in  that  direction.  If  the  initial  perturbation  was  instead  set  in  the  y  in¬ 
track  direction,  the  compensators  would  be  expected  to  estimate  and  control  them  better. 
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Settle  Time  (Orbits)  Control  (meters/second) 


Total  Compensator  Control  Over  One  Orbit,  Dy=10m 


Figure  47.  Compensator  Control,  Dy=10m,  Tuning  R’s 


Compensator  Settling  Times  Over  One  Orbit,  Dy=10m 
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Figure  48.  Compensator  Settling  Times,  Dy=10m,  Tuning  R’s 
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Compensator  Settling  Time  vs  Control  Over  One  Orbit,  Dy=10m 


Figure  49.  Compensator  Settling  Time  vs.  Control,  Dy=T0m,  Tuning  R’s 


The  AoTI  model  used  practically  the  same  control  as  AoTV,  but  it  allowed  longer 
settling  times  for  the  scenario.  Both  of  the  time-varying  models  decreased  their  settling 
times  for  this  model.  Table  9  below  compares  the  control  results  between  an  initial 
perturbation  in  the  x  radial  direction  versus  the  same  initial  perturbation  in  the  y  in-track 
direction.  The  unexpected  increase  in  control  for  this  scenario  is  explained  as  a  result  of 
cross-correlations  between  the  various  states  as  mentioned  in  the  previous  chapter.  For 
realistic  applications,  these  differences  would  be  statistically  studied  and  optimized  in 
diagonal  and  off-diagonal  weightings  of  filter  and  controller  tuning  matrices. 

Note  that  RMS  error  were  not  included  in  this  section  because  the  steady  state 
RMS  error  was  defined  earlier  as  a  zero  perturbation  calculation.  So  the  steady  state 
RMS  error  will  be  the  same  as  what  was  presented  in  Figure  43  of  Section  5.1.5.  The 
same  will  hold  true  for  next  two  scenarios,  Dz=10m  and  Dx=100m. 
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Table  9.  Dy  =  10  Scenario  Comparisons  for  R=lE-4 


Total  Control  Over  One  Orbit  (R=lE-4) 

Settling  Time  (R=lE-4) 

Control  (m/sec) 

%  over  DxlO  Scenario 

(%  orbit) 

%  over  DxlO 

AnTV 

0.0520 

10.3  % 

43.85 

-8.8  % 

AoTI 

0.0304 

4.4  % 

43.31 

5.5  % 

AoTV 

0.0304 

4.3  % 

41.21 

-0.5  % 

5.3  Initial  Perturbation  Dz  =  10  meters 

For  this  scenario,  the  initial  perturbation  of  10  meters  was  set  in  the  cross-track 
direction.  Again,  the  compensators  would  be  expected  to  estimate  and  control  them 
better  than  the  original  initial  radial  10  meter  perturbation. 


Total  Compensator  Control  Over  One  Orbit,  Dz=10m 


Figure  50.  Compensator  Control,  Dz=10m,  Tuning  R’s 
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Control  (meters/second)  Settle  Time  (Orbits) 


Compensator  Settling  Times  Over  One  Orbit,  Dz=10m 


Figure  51.  Compensator  Settling  Times,  Dz=10m,  Tuning  R’s 


Compensator  Settling  Time  vs  Control  Over  One  Orbit,  Dz=10m 


0.4  0.45  0.5  0.55 

Settle  Time  (Orbits) 


0.6 


0.65 


Figure  52.  Compensator  Settling  Time  vs.  Control,  Dz=10m,  Tuning  R’s 
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Similar  to  the  previous  scenario  but  this  time  more  obvious,  the  AoTI  and  AoTV 
had  practically  the  same  control  and  settling  time.  The  nodal  frame  time-varying  model 
actually  improved  settling  time  over  the  standard  Dx  =  10  m  scenario.  As  mentioned  for 
the  previous  scenario,  cross-correlations  between  states  are  the  culprits  responsible  for 
these  results.  Statistical  studies  of  these  behaviors  would  need  to  be  accomplished  and 
diagonal  and  off-diagonal  elements  of  the  tuning  matrices  would  be  altered. 


Table  10.  Dz  =  10  Scenario  Comparisons  for  R=lE-4 


Total  Control  Over  One  Orbit  (R=lE-4) 

Settling  1 

nine  (R=lE-4) 

Control  (m/sec) 

%  over  DxlO  Scenario 

%  orbit 

%  over  DxlO 

AnTV 

0.0504 

6.9  % 

44.22 

-8.1  % 

AoTI 

0.0333 

14.4  % 

41.41 

0.9  % 

AoTV 

0.0332 

14.0  % 

41.48 

0.2  % 

5.4  Initial  Perturbation  Dx  =  100  meters 

Although  the  initial  radial  perturbation  of  1 0  meters  would  be  a  realistic  scenario 
for  perturbations  that  could  occur  in  low  earth  orbit,  the  event  of  increased  solar  activity 
could  produce  even  larger  perturbations.  Increasing  the  initial  radial  perturbation  by  an 
order  of  magnitude  to  100  meters  would  not  produce  results  that  should  necessarily  be 
designed  to  for  stable  orbit  operations,  but  the  goal  for  the  compensators  should  be  that 
such  an  increased  perturbation  does  not  cause  them  to  diverge.  Orbital  reconfigurations 
as  mentioned  in  Chapter  I  would  require  the  controller  to  handle  perturbations  of  this 
magnitude.  We  would  expect  for  linear  systems  that  a  scaling  of  the  initial  perturbation 
should  produce  a  comparable  scaling  for  the  control.  This  scenario  will  test  the 
robustness  of  the  models. 
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Settle  Time  (Orbits)  Control  (meters/second) 


Total  Compensator  Control  Over  One  Orbit,  Dx=100m 


Figure  53.  Compensator  Control,  Dx=100m,  Tuning  R’s 


Compensator  Settling  Times  Over  One  Orbit,  Dx=100m 


Figure  54.  Compensator  Settling  Times,  Dx=100m,  Tuning  R’s 
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Compensator  Settling  Time  vs  Control  Over  One  Orbit,  Dx=100m 


Settle  Time  (Orbits) 


Figure  55.  Compensator  Settling  Time  vs.  Control,  Dx=100m,  Tuning  R’s 


While  the  magnitudes  for  control  jumped  an  order  of  magnitude,  the  general 
trends  seen  in  previous  scenarios  remained  the  same.  As  mentioned  at  the  beginning  of 
this  section,  for  linear  systems,  an  order  of  magnitude  increase  of  initial  perturbation 
should  produce  a  comparable  order  of  magnitude  increase  in  control.  Comparing  this 
scenario’s  results  with  the  standard  DxlO  scenario  (see  Table  1 1  below),  the  increase  in 
control  for  the  orbital  frame  models  over  the  standard  is  pretty  close  to  a  one-for-one, 
linear  increase.  A  10  times  increase  of  initial  perturbation  yielded  an  8.9  times  increase 
in  control.  But  the  nodal  frame  model  was  not  really  as  close.  For  AnTV  a  10  times 
increase  of  initial  perturbation  yielded  a  6  times  increase  in  control.  While  still  not 
besting  the  orbital  frame  models  in  control  usage,  AnTV  showed  that  the  varying  time- 
periodic  nature  for  its  dynamics  caused  less  deterioration  than  the  orbital  frame  models 
for  the  bigger  perturbation. 
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The  Aorbital  time-invariant  model  was  the  best  for  control,  and  it  also  provided 
the  best  settling  time.  Since  the  AoTI  and  AoTV  models  used  nearly  the  same  control, 
the  1%  better  settling  time  for  AoTI  would  make  it  a  better  option. 


Table  11.  Dx  =  100  Comparisons  for  R=lE-4 


Total  Control  Over  One  Orbit  (R=lE-4) 

Settling  Time  (R=lE-4) 

Control  (m/sec) 

%  over  DxlO 

Scenario 

(% 

orbit) 

%  over  DxlO 

Scenario 

AnTV 

0.3363 

613.5  % 

42.12 

76.5  % 

AoTI 

0.2890 

892.1  % 

41.15 

969.8  % 

AoTV 

0.2892 

892.1  % 

41.48 

629.1  % 

5.5  Initial  Perturbation  Dx  =  10  meters,  Increased  Step  Size 

For  the  purposes  of  this  thesis,  the  step  size  of  0.2338  seconds  (same  step  size 
used  in  Bordner’s  dynamics  model  [4])  has  been  ideal  for  ground  testing  applications,  but 
actual  space  computational  hardware  may  not  allow  or  want  this  measurement  frequency 
for  control  bums.  An  order  of  magnitude  increased  step  size  to  2.338  seconds  was  tested. 
While  sampling  theory  recommends  a  minimum  sample  time  2  to  3  times  the  frequency 
of  data  to  avoid  degradation  of  the  model,  this  scenario  tested  an  even  larger  increase. 
Again,  this  scenario  was  meant  not  to  produce  results  that  should  be  designed  to,  but 
more  as  a  test  to  see  if  the  higher  fidelity  models  would  diverge. 
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Avg  RMS  Error  (meters)  Control  (meters/second) 


Total  Compensator  Control  Over  One  Orbit,  Dx=10m,  Step  10 


Figure  56.  Compensator  Control,  Dx=10m,  Step  10,  Tuning  R’s 


Compensator  Average  RMS  Errors,  Dx=10m,  StepIO  (over  10  runs) 


Figure  57.  Compensator  RMS  Error,  Dx=0m,  Step  10  Tuning  R’s 
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Compensator  Settling  Times  Over  One  Orbit,  Dx=10m,  Step  10 


Figure  58.  Compensator  Settling  Time,  Dx=10m,  Step  10,  Tuning  R’s 

Since  this  scenario  incapacitated  the  time-invariant  model,  it  compounded 
evidence  of  the  robustness  of  both  the  AnTV  and  AoTV  models.  Even  with  its  hands  tied 
for  9/10  of  the  time,  AoTV  handled  the  scenario  with  only  6%  more  control  than  the 
standard  resulting  in  only  4%  better  settling  time.  The  AnTV  model  showed  marked 
control  usage  improvement  over  the  standard  step  size,  and  it  actually  improved  its  RMS 
error  by  over  31%  over  the  standard.  It  was  good  to  see  the  time-varying  models  handle 
this  as  well  as  they  did,  and  the  time-invariant  model  was  not  really  expected  to  perform 
well  in  this  type  of  environment  because  if  its  inherent  inflexibility.  Future  study  is 
recommended  to  look  at  the  step  size  limits  for  each  of  these  models.  For  AoTI  it  is 
recommended  to  evaluate  step  sizes  between  5  times  and  10  times  the  original  step  size. 
For  both  AnTV  and  AoTV,  evaluations  of  increased  step  size  above  the  10  times  increase 
are  recommended. 
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The  RMS  error  curves  in  Figure  57  show  that  for  the  steady  state  when  there  is  no 
initial  perturbation,  both  AoTV  and  AnTV  struggle  with  the  increased  time  step,  but 
AnTV  seems  to  handle  it  much  better. 


able  12.  Increased  Step  Size  Scenario  Comparisons,  DxMOm  for  R=lE-4 


Total  Control  (1  orbit) 

Settling  Time 

RMS  error  (Dx=0m) 

(m/sec) 

%  over  DxlO 

%  orbit 

%  over  DxlO 

(meters) 

%  over  DxlO 

AnTV 

0.0387 

-17.8  % 

30.31 

-37.0  %  * 

0.0231 

-31.5  % 

AoTI 

0.0381  * 

30.8  %  * 

>  100  * 

88.5  %  * 

0.0057 

385  % 

AoTV 

0.0310 

6.4  % 

40.57 

-3.6  % 

0.0552 

3430  % 

*  Diverging  system 


5.6  Initial  Perturbation  Dx  =  10  meters,  DGPS  Noise 

Earlier  scenarios  assumed  truth  model  measurement  noise  consistent  with 
absolute  GPS  position  readings.  Since  the  object  of  formation  flying  is  to  have  more  than 
one  satellite  flying  in  formation  with  a  number  of  other  cooperating  satellites,  the  use  of 
DGPS  should  be  a  viable  option.  According  to  Bordner’s  thesis,  the  DGPS  signals 
produce  accuracies  on  the  order  of  centimeters  [4].  For  this  scenario,  the  measurement 
noise  in  the  truth  model  was  changed  from  standard  deviation  of  Rx  =  4  meters  to  Rx  = 
0.04  meters  and  the  standard  deviations  for  Ry  and  Rz  were  similarly  reduced  two  orders 
of  magnitude  from  2  meters  to  0.02  meters.  Now,  the  reader  should  be  made  aware  that 
this  is  not  a  “true”  DGPS  scenario.  This  scenario  utilizes  advertised  capabilities  of 
DGPS  to  create  a  lower  noise  environment.  To  fully  evaluate  an  actual  DGPS  scenario, 
Bordner’s  thesis  provides  a  good  description  of  how  to  evaluate  this  more  complex 
problem  which  requires  multiple  cooperating  satellites  [4]. 
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Control  (meters/second) 


Total  Compensator  Control  Over  One  Orbit,  Dx=10m,  DGPS 


Figure  59.  Compensator  Control,  DGPS  Noise,  Tuning  R’s 
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Compensator  Settling  Times  Over  One  Orbit,  Dx=10m,  DGPS 


Figure  61.  Compensator  Settling  Time,  DGPS  Noise,  Tuning  R’s 


Compensator  Settling  Time  vs  Control  Over  One  Orbit,  Dx=10m,  DGPS 


Figure  62.  Compensator  Settling  Time  vs  Control,  DGPS  Noise,  Tuning  R’s 
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This  scenario  produced  results  consistent  with  expectations.  Note  that  RMS  error 
results  (Figure  60)  were  later  statistically  validated  over  100  Monte  Carlo  runs.  Reduced 
measurement  noise  for  the  system  reduced  both  the  control  and  RMS  error  for  all  models. 
But  the  nodal  frame  model  showed  the  most  improvement,  30%  better  control  and  98% 
better  RMS  error  as  compared  to  our  standard  noise  model.  Both  orbital  frame  models 
both  showed  modest  improvement  with  AoTV  again  having  the  best  RMS  error.  It  is  not 
surprising  that  AnTV,  with  its  previously-shown  strong  correlation  to  noise,  substantially 
reduced  its  RMS  error.  As  mentioned  at  the  beginning  of  this  section,  all  these  results  are 
for  a  simple  reduced-noise  scenario.  In  order  to  plan  for  operational  use  of  DGPS,  a  full- 
up  DGPS  scenario  using  methods  described  in  Bordner’s  thesis  would  be  needed  [4]. 


Table  13.  DGPS  Noise  Scenario  Comparisons,  Dx=10m  for  R=lE-4 


Total  Control  (1  orbit) 

Settling  Time 

RMS  error  (Dx=0m) 

(m/sec) 

%  over  DxlO 

%  orbit 

%  over  DxlO 

(meters) 

%  over  DxlO 

AnTV 

0.0329 

-30.2  % 

41.53 

-13.7% 

3.6E-4 

-97.9  % 

AoTI 

0.0288 

-1.2% 

41.26 

0.6  % 

2.9E-4 

-75.3  % 

AoTV 

0.0289 

-0.8  % 

41.50 

0.2  % 

1.7E-4 

-88.2  % 

5.7  Initial  Perturbation  Dx  =  10  meters,  Increased  Noise 

Earlier  scenarios  assumed  truth  model  measurement  noise  consistent  with 
absolute  GPS  position  readings.  Given  a  scenario  where  GPS  position  data  was  not 
available,  the  position  data  would  need  to  be  estimated  using  other  methods.  This 
scenario  (not  necessarily  the  worst  case  scenario)  would  increase  the  standard  deviation 
for  the  measurement  noise  by  two  orders  of  magnitude.  Similar  to  other  scenarios 
covered  in  previous  sections,  this  scenario  would  hopefully  never  be  needed  for  the 
actual  satellite  missions,  but  it  was  intended  as  a  test  of  the  fidelity  of  the  models. 
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Avg  RMS  Error  (meters)  £1  Control  (meters/second) 


Total  Compensator  Control  Over  One  Orbit,  Dx=10m,  lOOxNoise 


gure  63.  Compensator  Control,  Dx=10m,  lOOxNoise,  Tuning  R’s 


Compensator  Average  RMS  Errors  Over  One  Orbit,  Dx=0m,  lOOxNoise 


Figure  64.  Compensator  RMS  Error,  Dx=10m,  lOOxNoise,  Tuning  R’s 
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Compensator  Settling  Times  Over  One  Orbit,  Dx=10m,  lOOxNoise 


Figure  65.  Compensator  Settling  Time,  Dx=10m,  lOOxNoise,  Tuning  R’s 


Compensator  Settling  Time  vs  Control  Over  One  Orbit,  Dx=10m,  lOOxNoise 
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Figure  66.  Compensator  Settling  Time  vs.  Control,  Dx=10m,  lOOxNoise,  Tuning  R’s 
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Although  this  scenario  killed  the  nodal  model,  both  the  time-varying  and  time- 
invariant  models  in  the  orbital  frame  produced  acceptable  results  with  increased  control 
and  settling  time  as  compared  to  the  standard  (see  Table  14  below).  This  is  nice  to  see 
that  AoTI,  with  its  inherent  inadaptability,  was  able  to  handle  this  worse  case  scenario. 
AoTI’s  RMS  error  was  not  bad  with  only  4900%  more  error  than  the  standard  scenario, 
showing  that,  given  an  unperturbed  system,  the  time-invariant  model  still  worked  well. 

This  scenario  was  disappointing  for  the  AnTV  model,  because  it  had  perfonned 
fairly  well  over  a  wide  range  of  scenarios,  but  it  could  not  adapt  to  the  increased  noise, 
and  even  the  zero  perturbation  run  diverged.  These  past  two  scenarios  showed  a  strong 
correlation  between  measurement  noise  and  nodal-frame  model  perfonnance.  While  a 
satellite  system  would  not  necessarily  need  to  operate  over  such  a  wide  range  of 
measurement  noise  magnitudes,  this  performance  of  the  nodal  frame  model  definitely 
ruled  it  out  for  satellites  with  navigation  systems  that  produce  data  much  worse  than 
absolute  GPS  position  data.  On  the  other  hand,  this  was  good  news  for  the  AoTV 
because,  throughout  all  scenarios  tested,  AoTV  proved  its  fidelity  and  robustness  in  the 
face  of  many  challenges.  This  final  test  for  AoTV  produced  a  better  settling  time  than 
the  standard  and  the  smallest  increase  in  RMS  error  compared  with  the  other  two  models 
due  to  the  one  hundred-fold  increase  of  noise. 


Table  14.  Increased  Noise  Scenario  Comparisons,  Dx=10m  for  R=lE-4 


Total  Control  (1  orbit) 

Settling  Time 

RMS  error  (Dx=0m) 

(m/sec) 

%  over  DxlO 

%  orbit 

%  over  DxlO 

(meters) 

%  over  DxlO 

AnTV 

2.2403  * 

4653  %  * 

77.3  * 

107.9  %  * 

3.4412* 

19528  %  * 

AoTI 

0.0669 

129.6  % 

39.9 

88.5  % 

0.1650 

4900  % 

AoTV 

0.0455 

56.2  % 

44.6 

-3.6  % 

0.0324 

2080  % 

*  Diverging  system 
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VI.  Conclusions  and  Recommendations 


Several  conclusions  can  be  made  from  the  data  produced  by  this  thesis.  The 
earliest  lesson  learned  in  this  thesis  was:  do  not  try  to  solve  the  periodic  Riccati  equation 
using  units  of  the  earth’s  radius  expecting  to  get  centimeter-level  accuracy.  Numeric 
difficulties  arise  from  trying  to  invert  matrices  with  such  large  (1E4)  and  small  (IE- 17) 
elements.  Integrations  and  matrix  inversions  become  questionable  as  to  their  validity.  A 
possible  unit  conversion  could  possibly  cure  this  ill.  The  limits  of  numerical  precision 
for  computers  ranges  in  IE- 16  order  of  magnitude,  so  a  conversion  to  bigger  units  may 
avoid  crossing  the  minimum  numerical  precision  boundary. 

From  the  more  fruitful  SDRE  analysis,  data  shows  how  an  idealized  system  like 
the  full  state  feedback  regulator  may  not  always  predict  how  a  system  will  act  in  a  more 
realistic  noisy  environment.  Even  though  it  was  not  expected  to  perform  well,  the  AnTI 
model  looked  perfectly  adequate  for  the  regulator.  It  actually  produced  the  best  RMS 
error  compared  with  the  other  models  at  R=lE-4.  However,  once  you  took  away  its 
perfect  state  knowledge  and  only  gave  it  noisy  measurements,  it  provided  inadequate 
control  and  diverged.  This  model  diverged  for  all  output  feedback  compensator  scenarios 
including  the  scenario  with  measurement  noise  reduced  by  two  orders  of  magnitude. 

Even  in  a  low  noise  environment,  the  AnTI  model  diverged.  These  results  confirm  the 
limits  of  the  separation  principle  as  mentioned  in  Section  3.5.  The  robustness  of  the 
linear  quadratic  regulator  is  not  guaranteed  when  combined  with  the  Kalman  filter  to 
form  the  output  feedback  compensator.  As  mentioned  in  Chapters  4  and  5,  this  model 
was  not  expected  to  work  well,  but  this  is  a  lesson  on  how  results  from  idealistic  systems 
might  lead  you  down  the  wrong  path. 
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First  noted  for  the  full  state  feedback  regulator,  and  later  confirmed  with  most  of 
the  output  feedback  compensator  scenarios,  the  AoTI  and  AoTV  models  produced  very 
similar  results  for  most  of  the  scenarios.  Many  scenarios  produced  differences  in  both 
control  usage  and  settling  time  of  less  than  one  percent  difference  between  the  two.  This 
result  could  justify  the  use  of  linear  time  invariant  (LTI)  control  whose  validity  was 
questioned  in  Section  3.3.  Although  the  time-varying  control  produced  superior  results, 
the  similarity  in  results  combined  with  the  computational  capability  required  for  that 
method  may  drive  designers  to  favor  the  less  demanding  LTI  method.  But  the  failure 
scenarios  must  also  be  considered.  The  AoTI  model  diverged  in  the  face  of  reduced  step 
size.  This  failure  shows  the  inherent  inflexibility  of  the  LTI  method  in  compensating  for 
big  changes  over  large  intervals.  If  mission  requirements  require  control  at  a  time  step 
exceeding  one  second,  this  method  would  not  be  the  model  of  choice. 

The  nodal  frame  time-periodic  model  AnTV  showed  that  it  was  also  robust. 

AnTV  only  diverged  in  the  last  scenario,  but  it  consistently  produced  inferior  results  as 
compared  with  AoTV.  The  discrepancy  was  not  always  large,  but  AoTV  consistently 
produced  settling  times  and  control  4-7%  better  for  nominal  scenarios.  Another  factor  to 
consider  is  the  correlation  between  AnTV’s  performance  and  the  amount  of  measurement 
noise.  In  the  regulator  scenario  where  there  was  no  measurement  noise,  AnTV  actually 
produced  similar  RMS  error  to  AoTV.  For  the  DGPS  scenario  where  noise  was  reduced 
by  a  factor  of  100,  AnTV  was  only  5%  worse  in  RMS  error  than  AoTV.  Once  you  added 
more  realistic  noise,  AnTV  suffered.  The  one  hundred-fold  increase  in  noise  totally 
incapacitated  the  AnTV  model. 

The  correlation  between  the  periodic  nature  of  the  nodal  frame  A  matrix 
components  and  the  response  of  the  models  to  measurement  noise  was  speculated  in 
Section  5.1.5.  Basically  as  the  A  matrix  components  propagated  from  their  initial  values, 
the  measurement  noise  dominated  the  system  causing  the  control  to  vary  more  (compare 
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Figure  6  from  Section  2.3.1  with  Figure  35  from  Section  5.1.2).  As  time  approached  20- 
40%  of  the  orbit,  where  most  of  the  A  matrix  elements  periodically  increased  their 
values,  the  dynamics  started  to  dominate  and  the  system  was  controlled  with  less 
variability.  But  nearing  the  50%  mark  of  the  orbit,  most  of  the  non-zero  A  matrix 
components  again  approached  zero.  These  smaller  dynamics  components  allowed  the 
noise  again  to  dominate  giving  us  the  noisy  control. 

Correlation  between  noise  and  the  periodic  nature  of  the  dynamics  was  also  seen 
for  the  orbital  frame  time -varying  model,  but  it  was  to  a  much  lesser  degree  for  most 
scenarios.  For  the  scenario  where  the  measurement  noise  was  increased  by  two  orders  of 
magnitude,  a  similar  but  more  pronounced  trend  occurred  for  AoTV  (see  Figure  103  in 
Appendix  B).  Again  the  noise  dominates  the  system  at  the  initial  time,  50%  orbit,  and 
final  time  of  the  orbit.  Although  many  of  the  A  matrix  components  are  nearly  constant, 
some  components  have  the  same  periodic  nature  (just  a  much  smaller  magnitude)  as  we 
saw  in  the  nodal  frame.  As  these  components  go  to  zero,  the  noise  dominates.  Once  they 
periodically  increase  enough  (near  10%  or  65%  of  the  orbit),  the  dynamics  overshadows 
the  noise,  and  the  variable  control  is  noticeably  reduced.  Comparing  Figure  103  with 
Figure  105,  you  can  really  see  the  difference  between  time -varying  and  time-invariant 
control  in  this  more  noisy  environment.  For  AoTV  the  reduced  noise  areas  of  the  orbit 
allowed  the  model  to  control  in  a  more  productive,  less  noisy  way.  This  improvement 
gave  AoTV  50%  less  control  over  the  orbit  than  its  time-invariant  counterpart  AoTI. 

The  AoTV  model  showed  that  time  periodic  control  in  the  orbital  frame  produced 
not  only  the  best  results  for  most  scenarios,  but  it  was  the  only  one  never  to  diverge. 
AoTV  produced  sub-centimeter  level  accuracy  for  most  scenarios  using  less  than  3 
cm/sec  of  control  over  an  orbit.  If  the  scenario  was  carried  out  longer  than  one  orbit,  to 
cover  a  whole  day,  the  control  results  were  still  less  than  3.5  cm/sec  of  control.  If  this 
scenario  would  be  typical  for  every  day  of  the  year,  the  annual  control  usage  would  be 


81 


12.77  m/sec.  This  should  translate  to  become  an  acceptable  fuel  budget  for  our 
microsats.  These  results  are  more  impressive  when  you  note  that  the  noise  level  of  the 
system  was  30-50  times  worse  than  our  results.  In  other  words,  we  gave  the  system  noisy 
imperfect  data,  and  it  improved  that  data  by  3000-5000%. 

AoTV’s  higher  fidelity  would  be  a  good  selling  point  for  the  satellite  formation 
flying  mission,  because  the  mission  requires  inherent  flexibility.  Consistent  results 
through  various  scenarios  would  make  this  model  the  primary  contender  for  future  study. 
Before  this  model  could  be  considered  flight-worthy,  though,  much  more  tuning  and 
Monte  Carlo  simulation  runs  would  need  to  be  run  to  statistically  validate  its  capabilities. 
Further  tuning  would  need  to  look  at  the  position  state  cross-correlations  witnessed  by 
comparing  the  results  of  Dx,  Dy,  and  Dz  deflections.  Diagonal  and  off-diagonal 
weightings  for  the  filter  and  controller  tuning  matrices  would  have  to  be  considered, 
tested,  and  validated.  Further  tests  should  also  try  other  orbits  with  different  altitudes, 
eccentricities,  and  inclinations. 

This  thesis  focused  on  a  single  leader  satellite  and  the  control  techniques  in  order 
to  keep  the  satellite  set  in  the  desired  periodic  orbit.  For  a  more  populous  satellite 
cluster,  the  same  techniques  would  be  repeated  with  new  initial  conditions.  The  initial 
position  of  the  additional  satellite  would  be  chosen  to  provide  the  desired  cluster 
configuration.  As  the  satellite  propagates  through  its  orbit,  the  additional  satellite  is 
controlled  to  follow  its  own  prescribed  orbit  relative  to  the  leader.  Cooperative 
techniques  would  need  to  be  added  to  ensure  the  relative  spacing  and  positioning 
between  cluster  members  remains  adequate  for  the  desired  mission  at  all  times  during 
propagation  through  the  full  orbit. 
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Appendix  A.  Linear  Time  Periodic  Method  (LTP)  Troubleshooting 


As  mentioned  in  Section  3.4,  the  goal  of  the  Linear  Time  Periodic  (LTP)  method 
was  to  solve  the  time  periodic  Riccati  differential  equation.  Bittanti  proposed  the  method 
followed  in  this  thesis  with  many  checks  allowing  validation  of  the  results  [3]. 


A.l  Symplectic  Checks 


Following  the  method  detailed  in  Section  3.4,  checks  on  the  symplectic  nature  of 
certain  steps  were  used  to  help  validate  or  invalidate  the  results.  The  periodic  Riccati 
equation  we  were  trying  to  solve  is  Equation  36.  The  Hamiltonian  for  that  Riccati 
equation  is  given  by  Equation  37  (for  convenience  to  the  reader,  repeated  here  below) 


H(t)  = 


A(t)T 

-C(t)T  V2'1  C(t) 


-B(t)VlB(t)T 
—  A(t) 


(37) 


consisting  of  matrices  relating  to  the  state-space  relations  in  Equation  3 1 .  The 
first  check  was  a  symplectic  check  (Equation  38)  of  the  Hamiltonian.  Following  the 
periodic  nature  of  its  elements,  the  Hamiltonian  was  also  time  periodic,  so  the  check  was 
done  for  all  Hamiltonians  over  the  entire  orbit.  All  Hamiltonian  checks  were  successful, 
producing  zero  residuals. 

The  next  symplectic  check  was  not  successful.  The  Hamiltonian  transition  matrix 
<LH(t,t0)  (Equation  41)  should  follow  the  symplectic  relationship  in  Equation  40.  This  is  a 
critical  step  that  leads  directly  to  the  periodic  Riccati  equation  solution  (Equation  42),  so 
with  its  results  invalidated,  a  new  approach  would  have  to  be  taken  to  solve  our  problem. 
But  many  possible  trouble  areas  were  tested  first  to  see  where  our  calculations  broke 
down. 
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A.2  Integration  Methods 


Since  the  Hamiltonian  was  validated,  it  made  sense  to  look  into  the  integration 
method  used  to  calculate  <hH(t,t0),  which  was  the  integration  of 

4>H(t,to)  =  H-<EH(t,t0)  (51) 

As  mentioned  in  Section  2.2.1,  the  system  used  the  haming  integration  subroutine, 
an  ordinary  differential  equations  integrator  using  a  fourth  order  predictor-corrector 
algorithm.  Haming  collects  the  last  four  values  for  the  state  vector  and  extrapolates  them 
to  predict  the  next  value.  It  then  corrects  the  extrapolated  prediction  to  find  a  new  value 
for  the  state  vector.  Alternative  double  precision  canned  integration  methods  in  fortran’s 
IMSL  library  include  both  a  Runge-Kutta-Vemer  sixth-order  integration  routine  and  a 
Burlish-Stoer  integration  routine.  Both  of  these  routines  refused  to  handle  the  integration 
of  the  Hamiltonian  matrix,  claiming  the  system  was  “too  stiff’.  Haming  proved  to  be  the 
only  integrator  robust  enough  to  even  take  on  the  integration. 

Another  method  to  compute  <f>H(t,to),  using  Peano-Baker  series  expansion  [6],  was 
the  last  attempt 

t  t  C?! 

<J>H(t,t0)  =  I  + jHdaj  +  Jh- jHda.do,  +•••  (52) 

‘o  to  t0 

Again,  even  for  the  third  order  expansion,  the  computed  <f>H(t,t0)  still  failed  the 
symplectic  test  in  Equation  40. 

Of  course,  the  big  question  is  why  did  these  methods  not  work?  The  fortran 
routines  provide  a  good  hint  by  telling  us  the  system  is  “too  stiff’.  Trying  to  integrate  a 
12  x  12  dimensioned  matrix  with  elements  on  the  order  of  1E4  and  IE- 17  is  a  difficult 
task  especially  considering  the  limits  of  numerical  precision  for  computers  ranges  around 
IE- 16  order  of  magnitude.  Some  elements  of  the  Hamiltonian  are  already  pushing  the 
lower  limits  of  accuracy,  so  the  inaccurate  results  we  received  should  not  be  too  big  of  a 
surprise. 
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Appendix  B.  Selected  Data  Runs  for  Output  Feedback  Compensator 


B.l  Filter  Tuning  N1  and  N2 

After  tuning  the  R  controller  tuning  values  for  the  output  feedback  compensator, 
the  next  step  was  to  tune  the  N2  filter  tuning  values.  Similar  to  the  tuning  process  for  the 
controller’s  R  tuning  values,  the  object  is  to  minimize  both  control  and  RMS  error.  Also, 
as  before,  this  is  a  process  of  trade-offs.  Control  cost  increased  as  settling  time  improved 
with  lower  orders  of  magnitude  for  N2.  Steady  state  RMS  error  prefers  a  higher  N2 
value,  but  these  higher  N2’s  drive  the  system  divergent  when  a  perturbation  is  added. 
Further  tuning  of  the  N 1  values  (similar  to  the  Q  tuning)  was  not  as  big  a  factor,  due  to 
the  overbearing  nature  of  measurement  noise  over  dynamics  noise  for  our  system.  While 
N 1  tuning  values  were  also  quickly  looked  at,  results  were  similar  to  the  N2  results,  so 
the  end  result  left  us  with  our  original  tuning  values  N2=N1=1E0. 


Total  Compensator  Control  Over  One  Orbit,  R=1  E-4,  Dx=10m 


Figure  67.  Compensator  Control,  Dx=10m,  Tuning  N2’s 
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Settle  Time  (Orbits’ 


Compensator  Settling  Times  Over  One  Orbit,  Dx=10m 


Figure  68.  Compensator  Settling  Times,  Dx=T0m,  Tuning  N2’s 


Compensator  Average  RMS  Over  One  Orbit,  Dx=10m 


Figure  69.  Compensator  RMS  Error,  Dx=10m,  Tuning  N2’s 
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Although  the  trends  are  not  as  easy  to  distinguish,  each  model  has  its  own  filter 
tuning  value  N2  which  minimizes  its  control  or  settling  time.  For  each  orbital  frame 
model,  there  was  a  tuning  value  N2=lE-2  which  improved  settling  time  over  the  best  for 
our  initial  tuning  N2=1E0.  But  this  value  was  different  for  the  nodal  frame  model  with 
AnTV  favoring  N2=1E0.  For  an  example  of  AoTI  at  N2=lE-2,  a  0.2%  smaller  settling 
time  will  cost  4%  more  control.  While  control  could  possibly  be  improved  with  higher 
N2  values,  there  were  barriers  of  higher  N2’s  where  the  system  would  diverge  when  the 
filter  stopped  believing  noisy  measurement  data  and  only  believed  the  dynamics  model. 

The  steady  state  RMS  error  was  different  than  the  other  scenarios  because  it  has 
no  initial  perturbation.  With  no  perturbation,  it  preferred  a  higher  N2  value  basically 
telling  the  filter  to  ignore  the  measurements  and  only  believe  the  dynamics  model.  This 
is  a  situation  astrodynamicists  warn  about  for  Kalman  filters  [26].  For  unperturbed 
orbital  models  where  the  dynamics  are  well  known,  Kalman  filters  have  a  tendancy  to 
“go  to  sleep,”  reducing  the  weighting  on  the  measurements  to  zero.  This  situation  is 
acceptable  in  an  unperturbed  environment,  but  if  the  system  encounters  a  perturbation 
such  as  the  scenarios  in  this  thesis,  the  system  needs  to  use  both  dynamics  model  data 
and  noisy  measurement  data.  Since  the  initial  tuning  value  of  N2=1E0  produced  an 
acceptable  compromise  between  RMS  error  and  control,  it  was  used  as  the  standard  for 
comparisons  in  the  sections  that  follow.  The  equation  for  calculating  RMS  error  is  given 
below: 


RMSerror  = 


X  desired)  (y  Y  desired )  +  Z  desired) 

V  3 


(52) 
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Tab 

le  15.  Compensator  Total  Control  Over  One  Orbit,  Dx=10m,  Tuning  ] 

N2’s 

(m/s) 

N=lE-6 

N  =lE-5 

N  =lE-4 

N  =lE-3 

N  =lE-2 

N=1E-1 

N=1 

N=1E1 

AnTV 

0.0777 

0.0530 

0.0478 

0.0472 

0.0472 

0.0471 

0.0471 

0.0471 

AoTI 

0.0416 

0.0348 

0.0317 

0.0314 

0.0303 

0.0295 

0.0291 

0.289* 

AoTV 

0.0734 

0.0419 

0.0318 

0.0315 

0.0303 

0.0295 

0.0291 

0.0290 

*  Diverging  system 


Ta 


ble  16.  Compensator  Settling  Time  Over  One  Orbit,  Dx=10m,  Tuning  N2’s 


%  orbit 

N=lE-6 

N=lE-5 

N  =lE-4 

N=lE-3 

N  =lE-2 

N=1E-1 

N  =1 

N=1E1 

AnTV 

46.09 

45.34 

46.17 

47.27 

47.86 

48.05 

41.47 

48.13 

AoTI 

41.50 

41.53 

41.62 

41.06 

40.94 

40.98 

41.03 

>100* 

AoTV 

43.06 

42.55 

41.79 

41.33 

41.27 

41.34 

41.40 

41.47 

*  Diverging  system 


Table  17.  Compensator  Average  RMS  Error  Over  One  Orbit,  Dx=0m,  Tuning  N2’s 


Error  (m) 

N=TE-6 

N  =lE-4 

N  =lE-2 

N  =1E-1 

N  =1E0 

N  =1E1 

N  =1E2 

AnTV 

0.0446 

0.0202 

0.0170 

0.0353 

0.0175 

0.0167 

0.0108 

AoTI 

0.0237 

0.0123 

0.0045 

0.0038 

0.0030 

0.0016 

0.0028 

AoTV 

0.0393 

0.0103 

0.0044 

0.0012 

0.0015 

0.0011 

0.0011 

B.2  Initial  Perturbation  Dy  =  10  m 

The  first  scenario  was  for  a  10  meters  perturbation  in  the  y  in-track  direction. 


B.2.1  Time-invariant  Anodal  (AnTV). 


Control  (meters/second) 


Compensator  Error  (x  -  xdes),  Dy=10m,  R=1  E-4,  AnTV 


Figure  70.  Compensator  Error  for  R=lE-4,  Dy=10m,  AnTV 
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Compensator  Control,  Dy=10m,  R=1  E-4,  AnTV  (detail) 
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Figure  71.  Compensator  Control  for  R=lE-4,  Dy=10m,  AnTV 
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Control  (meters/second) 


B.2.2  Time-varying  Aorbital  (AoTV) 


Compensator  Error  (x-  Kdes),  Dy=10m,  R=1E-4,  AoTV 


Figure  72.  Compensator  Error  for  R=lE-4,  Dy=10m,  AoTV 


Figure  73.  Compensator  Error  for  R=lE-4,  Dy=10m,  AoTV 


90 


Control  (meters/second) 


B.2.3  Time-invariant  Aorbital  (AoTI) 


Compensator  Error  (x  -  xdes),  Dy=10m,  R=1E-4,AoTI 


Figure  74.  Compensator  Error  for  R=lE-4,  Dy=10m,  AoTI 


Compensator  Control,  Dy=10m,  R=1E-4,  AoTI  (detail) 


Figure  75.  Compensator  Control  for  R=lE-4,  Dy=10m,  AoTI 
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B.2.4  Model  Comparisons. 


Table  18.  Compensator  Total  Control  Over  One  Orbit,  Dy=10m,  Tuning  R’s 


control 

R=lE-6 

R=lE-5 

R=lE-4 

R=lE-3 

R=lE-2 

R=1E-1 

R=1E0 

(m/sec) 

AnTV 

0.0976 

0.0680 

0.0520 

0.0429 

0.0357 

0.0295 

0.0255 

AoTI 

0.0345 

0.0322 

0.0304 

0.0279 

0.0238 

0.0184 

0.0154 

AoTV 

0.0344 

0.0321 

0.0304 

0.0279 

0.0238 

0.0182 

0.0153 

Table  19.  Compensator  Settling  Times,  Dy=10m,  Tuning  R’s 


(% 

R=lE-6 

R=lE-5 

R=lE-4 

R=lE-3 

R=lE-2 

R=1E-1 

R=1E0 

orbit) 

AnTV 

43.74 

43.77 

43.85 

44.09 

44.62 

44.46 

61.74 

AoTI 

43.15 

43.19 

43.31 

43.67 

44.66 

46.64 

52.88 

AoTV 

41.09 

41.12 

41.21 

41.48 

42.18 

43.62 

51.53 

B.3  Initial  Perturbation  Dz  =  10  m 

The  next  scenario  was  for  a  10  meter  perturbation  in  the  z  cross-track  direction. 


B.3.1  Time-varying  Anodal  (AnTV). 


92 


Control  (meters/second) 


Compensator  Error  (x  -  xdes),  Dz=10m,  R=1E-4,  AnTV 


Figure  76.  Compensator  Error  for  R=lE-4,  Dz=10m,  AnTV 


Compensator  Control,  Dz=10m,  R=1  E-4,  AnTV  (detail) 


Figure  77.  Compensator  Control  for  R=lE-4,  Dz=10m,  AnTV 
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Control  (meters/second) 


B.3.2  Time-varying  Anodal  (AoTV) 


Compensator  Error  (x  -  xdes),  Dz=10m,  R=1  E-4,  AoTV 


Figure  78.  Compensator  Error  for  R=lE-4,  Dz=10m,  AoTV 


Figure  79.  Compensator  Control  for  R=lE-4,  Dz=10m,  AoTV 
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Control  (meters/second) 


B.3.3  Time-invariant  Anodal  (AoTI). 


Compensator  Error  (x  -  xdes),  Dz=10m,  R=1E-4,  AoTI 


Figure  80.  Compensator  Error  for  R=lE-4,  Dz=10m,  AoTI 


Compensator  Control,  Dz=10m,  R=1E-4,  AoTI  (detail) 


Figure  81.  Compensator  Control  for  R=lE-4,  Dz=10m,  AoTI 
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B.3.4  Model  Comparisons. 


Table  20.  Compensator  Total  Control  Over  One  Orbit,  Dz=10m,  Tuning  R’s 


control 

R=lE-6 

R=lE-5 

R=lE-4 

R=lE-3 

R=lE-2 

R=1E-1 

R=1E0 

(m/sec) 

AnTV 

0.0952 

0.0658 

0.0504 

0.0416 

0.0340 

0.0274 

0.0245 

AoTI 

0.0385 

0.0355 

0.0333 

0.0305 

0.0259 

0.0196 

0.0145 

AoTV 

0.0382 

0.0353 

0.0332 

0.0304 

0.0258 

0.0195 

0.0144 

Table  21.  Compensator  Settling  Times,  Dz=10m,  Tuning  R’s 


(% 

R=lE-6 

R=lE-5 

R=lE-4 

R=lE-3 

R=lE-2 

R=1E-1 

R=1E0 

orbit) 

AnTV 

43.74 

43.77 

43.85 

44.09 

44.62 

44.46 

61.74 

AoTI 

41.29 

41.39 

41.48 

41.73 

42.18 

40.55 

63.28 

AoTV 

41.36 

41.39 

41.48 

41.73 

42.18 

40.55 

63.28 

B.4  Initial  Perturbation  Dx  =  100  m 

The  next  scenario  was  for  a  perturbation  of  100  meters  in  the  x  radial  direction. 


B.4.1  Time-varying  Anodal  (AnTV). 
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Control  (meters/second) 


Compensator  Error  (x  -  xdes),  Dx=100m,  R=1E-4,  AnTV 


Figure  82.  Compensator  Error  for  R=lE-4,  Dx=100m,  AnTV 


Compensator  Control,  Dx=100m,  R=1E-4,  AnTV  (detail) 


Figure  83.  Compensator  Control  for  R=lE-4,  Dx=100m,  AnTV 
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Control  (meters/second) 


B.4.2  Time-varying  Anodal  (AoTV) 


Compensator  Error  (k  -  xdes),  Dx=100m,  R=1E-4,  AoTV 


Figure  84.  Compensator  Error  for  R=lE-4,  Dx=100m,  AoTV 


Compensator  Control,  Dx=100m,  R=1E-4,  AoTV  (detail) 


Figure  85.  Compensator  Control  for  R=lE-4,  Dx=100m,  AoTV 
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Control  (meters/second) 


B.4.3  Time-invariant  Anodal  (AoTI) 


Compensator  Error  (x  -  xdes),  Dx=100m,  R=1E-4,  AoTI 


Figure  86.  Compensator  Error  for  R=lE-4,  Dx=100m,  AoTI 


Compensator  Control,  Dx=100m,  R=1E-4,  AoTI  (detail) 


Figure  87.  Compensator  Control  for  R=lE-4,  Dx=100m,  AoTI 
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B.4.4  Model  Comparisons. 


Table  22.  Compensator  Total  Control  Over  One  Orbit,  Dx=100m,  Tuning  R’s 


control 

R=lE-6 

R=lE-5 

R=lE-4 

R=lE-3 

R=lE-2 

R=1E-1 

R=1E0 

(m/sec) 

AnTV 

0.3781 

0.3473 

0.3363 

0.3375 

0.3568 

0.3944 

0.4497 

AoTI 

0.3033 

0.2908 

0.2890 

0.2933 

0.3035 

0.3147 

0.3230 

AoTV 

0.3035 

0.2910 

0.2892 

0.2935 

0.3036 

0.3148 

0.3229 

Table  23.  Compensator  Settling  Times,  Dx=100m,  Tuning  R’s 


(% 

R=lE-6 

R=lE-5 

R=lE-4 

R=lE-3 

R=lE-2 

R=1E-1 

R=1E0 

orbit) 

AnTV 

42.00 

42.03 

42.12 

42.39 

43.04 

42.46 

34.76 

AoTI 

41.04 

41.06 

41.15 

41.41 

42.08 

43.15 

52.32 

AoTV 

41.36 

41.39 

41.48 

41.76 

42.47 

43.52 

52.88 

B.5  Initial  Perturbation  Dx  =  10  m,  Increased  Step  Size 

The  next  scenario  tested  was  for  an  order  of  magnitude  increase  in  step  size  with 
the  original  initial  perturbation  of  10  meters  in  the  x  radial  direction.  The  step  size  was 
increased  from  0.2338  seconds  to  2.338  seconds. 

B.5.1  Time-varying  Anodal  (AnTV). 
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Control  (meters/second) 


Compensator  Error  (x  -  xdes),  Dx=10m,  Step  10,  R=1E-4,  AnTV 


Figure  88.  Compensator  Error  for  R=lE-4,  Dx=10m,  Step  10,  AnTV 


Compensator  Control,  Dx=10m,  Step  10,  R=1E-4  ,  AnTV  (detail) 


Time  (%  orbit) 


Figure  89.  Compensator  Control  for  R=lE-4,  Dx=10m,  Step  10,  AnTV 
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Control  (meters/second) 


B.5.2  Time-varying  Anodal  (AoTV). 


Compensator  Error  (x  -  xdes),  Dx=10m,  Step  10,  R=1E-4,  AoTV 


Figure  90.  Compensator  Error  for  R=lE-4,  Dx=10m,  Step  10,  AoTV 


Figure  91.  Compensator  Control  for  R=lE-4,  Dx=10m,  Step  10,  AoTV 


102 


Control  (meters/second) 


B.5.3  Time-invariant  Anodal  (AoTI). 


Compensator  Error  (x  -  xdes),  Dx=10m,  Step  10,  R=1E-4,  AoTI 


Figure  92.  Compensator  Error  for  R=lE-4,  Dx=10m,  Step  10,  AoTI 


Compensator  Control,  Dx=10m,  Step  10,  R=1E-4,  AoTI  (detail) 


Time  (%  orbit) 


Figure  93.  Compensator  Control  for  R=lE-4,  Dx=10m,  Step  10,  AoTI 
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B.5.4  Model  Comparisons. 


Table  24. 

Compensator  Total  Control  Over  One  Orbit,  Dx=] 

0m,  Step 

10 

control 

R=lE-6 

R=lE-5 

R=lE-4 

R=lE-3 

R=lE-2 

R=1E-1 

R=1E0 

(m/sec) 

AnTV 

0.0476 

0.0405 

0.0387 

0.0421 

0.0738 

0.210  * 

0.173  * 

AoTI 

0.040  * 

0.039  * 

0.038  * 

0.037  * 

0.033  * 

0.028  * 

0.020  * 

AoTV 

0.0326 

0.0304 

0.0310 

0.0354 

0.0612 

0.135  * 

0.075  * 

*  Diverging  system 


Table  25.  Compensator  Settling  Times,  Dx=T0m,  Step  10 


(% 

R=lE-6 

R=lE-5 

R=lE-4 

R=lE-3 

R=lE-2 

R=1E-1 

R=1E0 

orbit) 

AnTV 

30.05 

29.84 

29.16 

87.47 

>100 

>100  * 

>100  * 

AoTI 

>100  * 

>100  * 

>100  * 

>100  * 

>100  * 

>100  * 

>100  * 

AoTV 

41.34 

41.18 

40.57 

38.44 

85.74 

>100  * 

>100  * 

*  Diverging  system 


Table  26.  Compensator  Average  RMS  Error,  Dx=0m,  Step  10 


(meters) 

R=lE-6 

R=lE-5 

R=lE-4 

R=lE-3 

R=lE-2 

R=1E-1 

R=1E0 

AnTV 

0.0263 

0.0263 

0.0263 

0.0263 

0.0263 

0.0263 

0.0263 

AoTI 

0.0057 

0.0054 

0.0057 

0.0058 

0.0058 

0.0060 

0.0063 

AoTV 

0.0975 

0.0808 

0.0552 

0.1098 

0.0014 

0.00810 

0.0008 
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B.6  Initial  Perturbation  Dx  =  10  m,  DGPS  noise 


The  next  scenario  tested  was  for  the  original  initial  perturbation  of  10  meters  in 
the  x  radial  direction,  but  now  the  truth  model  measurement  noise  was  decreased  in 
standard  deviation  by  two  orders  of  magnitude.  This  was  meant  to  simulate  the  benefit  of 
using  DGPS  position  data  instead  of  absolute  GPS  data. 


B.6.1  Time-varying  Anodal  (AnTV). 


Compensator  Error  (x  -  xdes),  Dx=10m,  DGPS,  R=1E-4,  AnTV 


Time  (%  orbit) 


Figure  94.  Compensator  Error  for  R=lE-4,  Dx=10m,  DGPS  noise,  AnTV 
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3  Compensator  Control,  Dx=10m,  DGPS,  R=1  E-4,  AnTV  (detail) 


Figure  95.  Compensator  Control  for  R=lE-4,  Dx=10m,  DGPS  noise,  AnTV 


B.6.2  Time-varying  Anodal  (AoTV). 


Compensator  Error  (x-  xdes),  Dx=10m,  DGPS,  R=1E-4,  AoTV 


Figure  96.  Compensator  Error  for  R=lE-4,  Dx=10m,  DGPS  noise,  AoTV 
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Figure  97.  Compensator  Control  for  R=lE-4,  Dx=10m,  DGPS  noise,  AoTV 


B.6.3  Time-invariant  Anodal  (AoTI). 


Compensator  Error  (x  -  xdes),  Dx=10m,  DGPS,  R=1E-4,  AoTI 


Figure  98.  Compensator  Error  for  R=lE-4,  Dx=10m,  DGPS  noise,  AoTI 
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Figure  99.  Compensator  Control  for  R=lE-4,  Dx=10m,  DGPS  noise,  AoTI 


B.6.4  Model  Comparisons. 


Table  27.  Compensator  Total  Control  Over  One  Or 


oit,  Dx=10m,  DGPS  noise 


Control-m/sec 

R=lE-6 

R=lE-5 

R=lE-4 

R=lE-3 

R=lE-2 

R=1E-1 

R=1E0 

AnTV 

0.0344 

0.0331 

0.0329 

0.0335 

0.0356 

0.0394 

0.0449 

AoTI 

0.0302 

0.0290 

0.0288 

0.0292 

0.0302 

0.0315 

0.0323 

AoTV 

0.0304 

0.0291 

0.0289 

0.0293 

0.0304 

0.0315 

0.0323 

Table  28.  Compensator  Settlin 


g  Times,  Dx=10m,  DGPS  noise 


% 

orbit 

R=lE-6 

R=lE-5 

R=lE-4 

R=lE-3 

R=lE-2 

R=1E-1 

R=1E0 

AnTV 

41.41 

41.44 

41.53 

41.80 

42.44 

41.91 

34.58 

AoTI 

41.14 

41.40 

41.26 

41.52 

42.21 

43.30 

52.15 

AoTV 

41.38 

41.40 

41.50 

41.77 

42.48 

43.54 

52.87 
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Table  29.  Compensator  Average  RMS  Error,  Dx=0m,  DGPS  noise 


(meters) 

R=lE-6 

R=lE-5 

R=lE-4 

R=lE-3 

R=lE-2 

R=1E-1 

R=1E0 

AnTV 

5.3E-4 

2.4E-4 

3.6E-4 

8.3E-4 

4.3E-4 

3.9E-4 

4.3E-4 

AoTI 

3.3E-4 

3.5E-4 

2.9E-4 

3.5E-4 

4.4E-4 

5.6E-4 

4.6E-4 

AoTV 

2.5E-4 

8.9E-4 

1.8E-4 

6.5E-4 

1.2E-4 

2.2E-4 

4.7E-4 

B.7  Initial  Perturbation  Dx  =  10  m,  100  x  noise 

The  last  scenario  tested  was  for  the  original  initial  perturbation  of  10  meters  in  the 
x  radial  direction,  but  now  the  truth  model  measurement  noise  was  increased  in  standard 
deviation  by  two  orders  of  magnitude.  This  was  meant  to  simulate  the  detriment  of 
losing  the  ability  to  collect  absolute  GPS  data. 

B.7.1  Time-varying  Anodal  (AnTV). 
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Compensator  Error  (x  -  xdes),  Dx=10m,  lOOxNoise,  R=1E-4,  AnTV 


Figure  100.  Compensator  Error  for  R=lE-4,  Dx=10m,  100  x  noise,  AnTV 


Compensator  Control,  Dx=10m,  lOOxNoise,  R=1E-4,  AnTV 
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Figure  101.  Compensator  Control  for  R=lE-4,  Dx=10m,  100  x  noise,  AnTV 
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B.7.2  Time-varying  Anodal  (AoTV). 


Compensator  Error  (k  -  xdes),  Dx=10m,  lOOxNoise,  R=1E-4,  AoTV 


Figure  102.  Compensator  Error  for  R=lE-4,  Dx=10m,  100  x  noise,  AoTV 


Compensator  Control,  Dx=10m,  lOOxNoise,  R=1E-4,  AoTV  (detail) 


Figure  103.  Compensator  Control  for  R=lE-4,  Dx=10m,  100  x  noise,  AoTV 
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Control  (meters/second) 


B.7.3  Time-invariant  Anodal  (AoTI). 


Compensator  Error  (x-  xdes),  Dx=10m,  lOOxNoise,  R=1E-4,  AoTI 


Figure  104.  Compensator  Error  for  R=lE-4,  Dx=10m,  100  x  noise,  AoTI 


Compensator  Control,  Dx=10m,  lOOxNoise,  R=1E-4,  AoTI  (detail) 


Time  (%  orbit) 


Figure  105.  Compensator  Control  for  R=lE-4,  Dx=10m,  100  x  noise,  AoTI 
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B.7.4  Model  Comparisons. 


Table  30.  Compensator  Total  Control  Over  One  Orbit,  Dx=10m,  100  x  Noise 


control 

R=lE-6 

R=lE-5 

R=lE-4 

R=lE-3 

R=lE-2 

R=1E-1 

R=1E0 

(m/sec) 

AnTV 

6.983  * 

3.936  * 

2.240  * 

1.306  * 

0.793  * 

0.530  * 

0.426  * 

AoTI 

0.1728 

0.1038 

0.0669 

0.0488 

0.0399 

0.0356 

0.0334 

AoTV 

0.0901 

0.0607 

0.0455 

0.0379 

0.0340 

0.0331 

0.0327 

*  Diverging  system 


Table  31.  Compensator  Settling  Times,  Dx=10m,  100  x  Noise 


(% 

R=lE-6 

R=lE-5 

R=lE-4 

R=lE-3 

R=lE-2 

R=1E-1 

R=1E0 

orbit) 

AnTV 

>100  * 

>100  * 

>100  * 

>100  * 

>100  * 

>100  * 

>100  * 

AoTI 

77.26 

77.28 

77.34 

77.56 

78.34 

81.34 

93.96 

AoTV 

39.74 

39.78 

39.90 

40.28 

41.22 

42.56 

61.44 

*  Diverging  system 


r 

fable  32.  Compensator  Average 

RMS  Error,  Dx=0m,  1 

[00  x  Noise 

(% 

R=lE-6 

R=lE-5 

R=lE-4 

R=lE-3 

R=lE-2 

R=1E-1 

R=1E0 

orbit) 

AnTV 

1.5859  * 

4.2933  * 

3.4412  * 

4.7513  * 

2.5248  * 

3.8470  * 

2.1767  * 

AoTI 

0.0737 

0.0864 

0.1650 

0.1162 

0.1623 

0.2044 

0.2268 

AoTV 

0.1088 

0.1798 

0.0324 

0.1225 

0.3212 

0.1213 

0.0870 

*  Diverging  system 
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